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Abstract 


The  need  to  meet  the  demand  for  high  quality  digital  images,  with  comparatively  modest 
storage  requirements,  is  driving  the  development  of  new  image  compression  techniques. 
This  demand  has  spurred  new  techniques  based  on  time  to  frequency  spatial  transformation 
methods.  At  the  core  of  these  methods  are  a  family  of  transformations  built  on  basis  sets 
called  “wavelets.”  The  wavelet  transform  permits  an  image  to  be  represented  in  a 
substantially  reduced  space  by  transferring  the  energy  of  the  image  to  a  smaller  set  of 
coefficients.  Although  these  techniques  are  lossy  as  the  compression  ratio  rises,  very 
adequate  reconstructions  can  be  made  from  surprisingly  small  sets  of  coefficients.  This  work 
explores  the  transformation  process,  storage  of  the  representation  and  the  application  of  these 
techniques  to  24-bit  color  images.  A  working  color  image  compression  model  is  illustrated. 
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Chapter  1 


Introduction 


It  is  apparent  that  we  are  on  the  verge  of  another  revolution.  Not  unlike  the  personal 
computer  industry  of  only  five  or  ten  years  ago  influencing  the  way  that  work  is 
accomplished  in  the  office  environment  today,  we  are  now  facing  a  leap  in  the  way 
information  is  presented  to  the  consumer.  This  has  broad  implications  across  the  board. 
From  encyclopedias,  magazines  and  other  published  media,  to  interfaces  on  the  Internet,  to 
leisure  activities  such  as  games  and  video  taped  movies,  to  the  art  of  photography,  are  all 
being  influenced  by  the  ability  to  digitize,  store,  and  retrieve  images.  Images  are  the  most 
powerful  way  to  communicate  ideas  and  information  to  another  person,  not  only  by  crossing 
cultural  and  linguistic  barriers,  but  by  providing  a  common  frame  of  reference  for  the 
communication  of  ideas. 

During  the  past  ten  years  we  have  increased  our  ability,  at  the  consumer  level,  to  store 
and  access  information  by  1,000  folds.  At  the  same  time,  the  demand  on  our  capability  to 
store  and  access  this  information  has  exceeded  the  1,000  fold  gains.  New  demands  will 
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probably  double  or  triple  the  throughput  requirements  for  large  amounts  of  image  data  in  real 
time  applications  over  the  next  five  years. 

A  simple  example  of  the  kind  of  requirement  needed  in  the  near  future  is  digital  video. 
A  broadcast  quality  video  picture  might  have  512  X  512  pixels  and  24  bit  (RGB)  color 
resolution.  This  gives  us  a  digital  image  of  786,432  bytes  with  16.8  million  colors. 
Television  pictures  are  synced  at  30  frames  per  second  to  give  us  continuity  of  motion.  So 
to  digitize  one  second  of  the  video  frame  signals  requires  more  than  23.5  megabytes  of 
transmission  and  storage  capability.  One  hour’s  worth  of  digitized  video  would  require  84.9 
gigab5^es  of  storage  and  transmission  capability.  Most  personal  computers  today  have  500 
megabytes  of  storage  giving  a  capability  of  about  22  seconds  of  raw  video  storage.  The  CD- 
ROM  is  capable  of  storing  about  650  megabytes  of  data,  but  to  store  an  hour  of  video  images 
would  require  a  compression  ratio  of  over  1 30;  1 . 

Using  this  example  as  a  starting  point,  the  digital  representation  of  an  image  requires 
a  significant  number  of  bytes.  The  goal  of  image  compression  is  to  reduce  this  number  as 
much  as  possible  while  retaining  the  ability  to  restore  a  faithful  reconstmction  of  the  original 
image.  There  is  a  tradeoff  involved  in  the  exactness  of  the  reconstruction  and  level  of 
compression  desired.  No  matter  the  technique  involved  in  compression  there  at  some  point 
(generally  around  2: 1  compression  ratio)  begins  a  degradation  of  the  restored  image.  Then 
of  course,  a  compromise  between  the  acceptable  loss  of  fidelity  and  achievable  compression 


must  be  reached. 
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Every  image  capturing  system,  be  it  a  film  camera,  or  a  video  camera,  produces  an 
image  by  sampling  in  space,  quantizing  in  brightness,  and  transforming  into  a  color 
representation  the  analog  scenes.  A  photograph  does  this  by  modifying  the  emulsion  of  the 
film  (actually,  the  discrete  silver  halide  crystals)  and  a  video  camera  by  digitizing  into  pixels 
the  sampled  values.  The  sampling  step  size  is  usually  chosen  to  be  small  enough  to  take 
advantage  of  the  ability  of  the  human  visual  system  to  integrate  the  data  into  a  continuous 
image. 


Problem  Definition 

We  are  therefore  faced  with  a  problem  on  how  to  meet  the  demand  of  efficient  storage. 
Today’s  JPEG  standard  does  not  provide  optimal  performance  as  we  will  see.  JPEG  [32] 
stands  for  the  Joint  Photographic  Experts  Group.  It  is  referred  as  a  "joint"  group  because  this 
committee  is  sanctioned  by  the  CCITT  and  the  ISO,  two  prominent  international  standard 
groups.  JPEG  refers  both  to  the  committee  and  their  compression  standard  that  defines  a 
method  for  compressing  photographic  images.  Images  compressed  with  the  JPEG  algorithm 
undergo  a  "lossy"  compression.  The  amount  of  compression  can  be  varied,  with  a  resultant 
loss  or  gain  in  resolution.  JPEG  compression  can  achieve  impressive  compression  ratios 
compared  to  the  statistical  models,  reducing  the  storage  required  by  images  to  less  than  10% 
of  the  size  of  the  original  with  only  very  slight  loss  of  resolution.  By  sacrificing  more 
resolution  compression  can  reach  to  95%  or  more  using  JPEG. 
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The  problem  therefore  is  to  develop  a  compression  scheme  competitive  with  JPEG  that 
will  restore,  with  acceptable  loss,  an  image  while  maximizing  its  storage.  To  this  end  we 
want  to  develop  a  technique  that  will  produce  compression  ratios  of  20:1  or  better  for  color 
image  storage. 


Constraints  and  Assumptions 

Let  us  examine  the  properties  of  images  and  see  how  this  relates  to  image  compression: 
First,  data  in  an  image  is  not  random.  Adjacent  samples  have  similar  values  usually.  This 
correlation  or  “redundancy”  relates  to  the  statistical  properties  of  an  image  and  is  a  function 
of  resolution,  bit-depth,  image  noise  and  image  detail.  If  we  exploit  this  redundancy  then 
there  is  no  doubt  that  we  can  reduce  the  number  of  bits  required  for  the  original  image. 
Secondly,  there  is  a  component  of  “irrelevancy”  which  relates  to  the  observer’s 
comprehension  of  an  image  and  depends  on  image  noise,  detail,  and  viewing  conditions. 

We  want  a  method  that  is  flexible  enough  to  manage  images  of  varying  sizes  and 
complexity.  We  will  assume  that  these  images  convey  meaningful  information  from  the 
average  viewer’s  perspective  and  are  therefore  neither  white  noise  nor  periodic  in  nature,  but 
somewhere  in  between. 

First,  let  us  consider  the  basic  representation  of  data.  We  will  assume  a  message  to  be 
a  data  set  representing  either  text  or  an  image.  Although  we  are  concentrating  on  image 
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representation,  most  of  the  techniques  can  also  be  applied  to  text.  The  alphabet  of  the 
message  is  the  set  of  all  of  the  possible  symbols  that  may  appear  in  a  message  or  image.  For 
example,  when  compressing  ASCII  text  files,  the  alphabet  would  probably  consist  of  127 
characters,  0x00  through  0x7f.  An  image  in  8-bit  greyscale  format  would  have  an  alphabet 
of  256  symbols,  0x00  through  Oxff,  while  a  24-bit  RGB  color  image  has  an  alphabet  of 
nearly  16.8  million  symbols. 

Compression  schemes  normally  maintain  a  model,  which  is  a  set  of  accumulated 
statistics  describing  the  state  of  the  encoder.  For  example,  in  a  simple  compression  program, 
the  model  may  be  a  count  of  the  frequency  of  every  symbol  in  an  input  file 

What  is  the  minimum  number  of  bits  that  can  be  used  to  represent  a  message?  To  put 
this  in  context  let  us  examine  the  notion  of  entropy.  Entropy  is  the  measure  of  the  amount 
of  order  in  a  message,  a  small  value  means  there  is  a  great  deal  of  redundancy  and  a  large 
value  suggests  a  great  deal  of  disorder.  Entropy  therefore  is  a  measure  of  the  information 
content  of  a  message  with  respect  to  a  particular  model.  It  also  determines  the  minimum 
number  of  bits. per  symbol  needed,  on  average,  to  represent  a  message  generated  by  that 
model.  Different  coding  models  can  produce  different  probability  estimates.  For  a  model 
that  assigns  a  constant  probability  to  each  symbol  of  the  alphabet,  the  entropy  is  calculated 
as 

N 

E  =  -'LPi  loga  Pi 


(1) 
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where  /?,  is  the  probability  of  symbol  /.  Since  the  entropy  of  a  message  represents  the  lower 

bound  of  the  number  of  bits  required  to  losslessly  represent  the  message,  the  only  way  to 
exceed  this  limit  is  to  either  change  the  probability  generating  model  or  to  introduce  loss  into 
the  message. 


Chapter  2 


Related  Research 

Information  and  Models 

Even  before  the  advent  of  computers  there  were  attempts  to  minimize  the  alphabet 
used  to  represent  a  message.  Morse  code  is  an  early  example  of  a  variable  length  coding 
scheme  that  reduced  the  code  size  for  more  common  symbols.  Reducing  the  symbol  set  from 
a  fixed  representation  to  a  variable  representation  permits  information  to  be  encoded  more 
efficiently.  More  recently,  there  has  been  a  wealth  of  research  on  the  reduction  of  message 
representation  requirements  [1,4,13,12,14,15,16,22,26,34,37]. 

In  1948,  C.  E.  Shannon  [29]  and  R.  M.  Fano  independently  developed  a  coding 
technique  that  attempted  to  minimize  the  number  of  bits  used  to  encode  a  message  when  the 
probabilities  of  symbols  in  the  message  were  known.  Shortly  afterward,  D.  A.  Huffman  [15] 
developed  a  technique  that  superseded  Shannon-Fano  coding  and  produces  provably 
optimum  code  sets,  resulting  in  marginally  better  performance  than  Shannon-Fano  codes. 
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Huffman  code  is  less  than  optimal  when  the  probabilities  of  the  symbol  are  not  powers  of 
two. 


Another  variation  on  the  Huffman  coding  is  run-length  encoding,  where  the  count  of 
the  run  of  a  particular  symbol  is  included  in  coding  the  scheme.  If  the  symbol  is  optimized 
based  on  the  number  of  unique  runs  and  the  runs  are  relatively  long,  this  process  can  also 
reduce  the  number  of  symbols  required  for  coding  a  message.  This  is  particularly  well  suited 
for  binary  images  since  there  are  only  two  symbols  and  the  runs  can  be  very  long.  As  the 
number  of  symbols  increases  and  the  runs  become  very  short,  this  is  a  less  efficient 
technique. 

In  1987,  a  new  more  efficient  method  was  introduced  by  I.  H.  Witten,  R.  M.  Neal,  and 
J.  G.  Cleary  [34].  This  method  called  arithmetic  coding,  works  by  assigning  codes  to 
symbols  that  have  a  known  probability  distribution  and  takes  an  entire  message  and  encodes 
it  as  a  single  floating  point  number  less  than  1  and  greater  than  or  equal  to  0.  It  achieves  an 
average  code  length  arbitrarily  close  to  the  entropy  of  the  message.  Arithmetic  coding  can 
more  efficiently  encode  texts  or  images  by  eliminating  the  quantization  effects  of  other 
coding  techniques. 

These  encoding  techniques  are  lossless  and  only  reduce  the  space  required  by 
minimizing  the  representation  alphabet.  As  end  techniques,  they  can  also  be  fed  probabilities 
from  another  model.  These  are  more  sophisticated  models  that  convert  the  original  data  into 
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another  representation  that  improves  compression.  We  will  look  at  some  probability  models 
and  how  they  generate  the  statistics  that  are  fed  to  the  encoder. 

There  are  two  approaches  to  image  compression.  Techniques  based  on  exploitation 
of  predictive  means  of  an  image  are  called  “spatial  coding”  techniques.  The  first  group  of 
techniques  exploits  the  redundancy  in  the  data  and  attempt  to  predict  the  image  causally.  For 
example,  an  image  with  one  color  is  fully  predictable  once  the  first  pixel  is  known. 
Conversely,  a  white  noise  image  is  totally  unpredictable  and  every  pixel  has  to  be  stored  to 
reconstruct  the  image.  These  forms  of  compression  can  be  viewed  as  the  attempt  to 
transform  the  original  image  array  {u. .} ,  into  another  array  {v. .} ,  which  has  no  redundancy 
and  such  that  { u. . }  can  be  uniquely  determined  from  { v. .  j .  The  raw  data  rate  of  { v . . }  then 
determines  the  data  rate  of  {«,  ■}.  Generally,  this  type  of  process  results  in  compression  but 
with  an  accompanying  distortion  in  the  reconstructed  array  {«',  ,} . 

The  second  approach  achieves  compression  by  transforming  a  given  image  into  another 
array  such  that  the  maximum  information  is  packed  into  a  minimum  number  of  samples. 
These  are  based  on,  what  would  be  loosely  called,  a  frequency  domain  process.  This  is 
known  as  “transform  coding  ”  and  are  related  to  the  non-causal  representation  of  signals. 
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Spatial  Methods 

Lets  briefly  consider  the  predictive  coding  techniques.  The  principal  is  very  simple. 
Since  the  image  source  is  assumed  to  be  highly  correlated,  then  for  any  given  element  its 
neighbors  should  be  closely  related  in  value.  So,  we  can  use  the  preceding  values  of 
neighbors  from  the  same  line  or  previous  lines  to  predict  the  value  of  the  present  element. 
Given  our  expectation  of  a  highly  correlated  image,  the  difference  of  the  predicted  value  and 
the  actual  value  should  be  very  small.  The  only  error  introduced  is  in  quantizing  the 
difference  signal.  The  weakness  of  the  predictive  methods  is  that  they  are  very  sensitive  to 
variations  in  the  input  data  statistics  and  if  a  coding  error  occurs,  it  propagates  throughout 
the  remainder  of  the  message. 

Pulse-code  modulation  (PCM)  [4].  Standard  pulse  code  modulation  (PCM)  encoding 
is  a  common  technique  for  encoding  audio  data.  Telephone  conversations  and  audio  CDs 
both  use  conventional  PCM.  PCM  samples  a  waveform  at  uniform  steps  and  encodes  the 
level  of  the  waveform.  Acceptable  quality  images  can  be  obtained  from  compression  ratios 
of2.6.T. 

Differential  pulse-code  modulation  (DPCM)  [13].  DPCM  does  not  encode  the  level 
as  PCM  does,  it  instead  encodes  the  difference  from  the  last  sample.  Adaptive  DPCM 
(ADPCM)  [12]  takes  that  a  step  further,  and  modifies  the  coding  of  the  difference  depending 
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on  the  state  of  the  waveform.  This  technique  generates  ratios  about  2.5;  1,  but  adaptive 
versions  can  obtain  ratios  as  high  as  3.5: 1 . 

Interpolative  coding  [4,23].  Most  interpolative  coders  are  of  the  zero-order  or  first- 
order  models,  giving  compression  ratios  of  around  4:1.  Higher  order  polynomials  can  be 
used,  but  the  computational  complexity  generally  offsets  the  gain  in  the  compression  ratios. 

Bit-plane  coding  [4].  This  technique  uses  each  bit  of  a  byte  as  a  parallel  map  of  ones 
and  zeros,  these  are  then  run-length  encoded.  Run-length  encoding  encodes  the  symbol  and 
then  the  count  of  the  number  of  consecutive  repetitions  of  the  symbol.  Compression  ratios 
of  4: 1  can  be  obtained  without  using  any  visual  system  considerations. 

Prediction  by  Partial  Matching  (PPM)  [4].  This  is  a  technique  based  on  a  finite 
context  model.  PPM  makes  its  prediction  based  on  a  dynamic  Markov  model  of  various 
orders,  most  commonly  order- 1  or  order-2.  The  limiting  factor  is  the  trade  off  between 
additional  gains  made  in  compression  versus  the  space  requirement  of  this  model.  Each 
order  i  of  the  model  requires  M  possible  strings  to  be  considered,  where  N  is  the  number  of 
symbols  in  the  alphabet.  Compression  ratios  of  4: 1  can  be  achieved. 

Dictionary-based  methods  [4].  These  may  be  either  static  or  adaptive.  Macro 
substitution  methods  use  a  dictionary  to  compress  data.  A  string  of  symbols  is  encoded  as 
a  pointer  into  a  codebook  or  dictionary.  A  static  dictionary  will  compress  an  entire  stream 
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using  the  same  dictionary.  An  adaptive  method,  such  as  one  of  the  Ziv-Lempel  [37] 
schemes,  is  continually  modifying  its  dictionary.  Compression  ratios  as  high  as  8:1  can  be 
obtained. 


Transform  Methods 

It  is  at  this  point  we  will  leave  text  compression  behind,  since  text  compression 
requires  exact  reconstruction  and  the  following  methods  are  subject  to  finite  arithmetic 
precision  errors.  An  alternative  to  predictive  coding  is  transform  coding  where  the  operation 
is  a  successive  multiplication  of  elements  by  a  set  of  basis  vectors  {u„  U2,  ...,  u^,].  Basis 
vectors  in  this  case  allow  any  solution  to  expressed  uniquely  as  a  linear  combination  of  { Uj, 
U2, M„, }  in  V(R)  the  space  of  all  square  integrable  functions. 

Karhunen-Loeve  transform  (KLT):  This  is  the  best  linear  transformation  method. 
It  is  obtained  by  determining  the  basis  vector  set  which  diagonalizes  the  specific  data 
covariance  matrix,  and  it  results  in  a  transform  domain  covariance  matrix  of  uncorrelated 
components.  This  packs  the  maximum  average  energy  into  a  minimal  set  of  samples.  Since 
there  is  no  fast  algorithm  associated  with  it,  it  is  not  used  because  of  the  computational  load 
required. 
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Fast  transformations:  These  transforms  take  a  set  of  data  in  the  spatial  domain  and 
convert  them  into  an  identical  representation  in  the  frequency  domain.  An  important 
difference  between  these  transforms  and  the  KLT  is  that  these  transforms  do  not  depend  on 
the  input  image  statistics.  The  three  important  features  of  a  suitable  transform  are  its 
compressional  efficiency,  which  relates  to  concentrating  the  energy  at  low  frequency,  ease 
of  computation,  and  minimum  mean  square  error.  One  of  the  key  features  of  the  fast 
transforms  is  that  they  can  be  computed  with  N  log2  N  or  fewer  operations  as  compared  to 
in  the  standard  Fourier  transform. 

Discrete  Fourier  transform  (DFT):  The  Fourier  transform  converts  the  signal 
from  the  time  domain  to  the  frequency  domain.  This  new  domain  has  basis  functions  that 
are  cosines  and  sines.  The  drawbacks  to  the  Fourier  transform  is  that  it  poorly  localizes 
components  in  time  and  the  transform  is  represented  by  both  a  real  and  imaginary 
component.  It  is  described  by: 

X(k)  =  ^  k  =  0,...,A^-1  (2) 

n=0 


Discrete  cosine  transform  (DCT)[2)1\.  This  has  many  of  the  same  properties 
of  the  KLT  and  it  can  be  thought  of  as  consisting  of  the  cosine  portion  (the  real  part)  of  the 
DFT.  The  DCT  is  used  as  the  core  of  the  Joint  Photographic  Experts  Group  (JPEG) 
standard.  It  is  represented  by: 
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N-l 


X(k)  = 


cos 


k'sz{ln  + 1 ) 


IN 


,  k  - 


(3) 


Discrete  wavelet  transform  (DWr}[36\:  In  the  wavelet  domain  the  basis 
functions  are  more  complicated  than  the  sine  and  cosine  functions  of  the  DFT.  These  basis 
functions  are  called  ’’mother  functions”  or  “wavelets.”  Unlike  the  sines  and  cosines  which 
define  the  Fourier  transform,  there  is  not  a  single  unique  set  of  wavelets.  It  is  this  last  class 
of  fast  transform  methods  that  we  will  use  as  the  core  of  our  compression  method. 


Chapter  3 


Preliminaries 


Transforms 


In  order  for  us  to  understand  how  the  wavelet  transforms  work,  we  need  to  discuss  how 
transforms  work  in  general.  The  definition  of  a  Mh  order  linear  transform  of  a  one 
dimensional  sequence  x(n);  n  =  0,  I,  is  given  by: 

N-l 

yik)  =  Y,  x{n)fik,n)  for  k  =  0,l,-,N-l  (4) 

n=0 


'Where  J{k,n)  is  a  forward  transformation  kernel,  and  y(^)  are  the  transform  coefficients.  The 
inverse  transform  that  recovers  the  input  sequence  is 

N-\ 

x(n)  =  Y,y(k)g{k,n)  for  n  =0, 1,...,A^-1  (5) 

k=0 

where  g(k,n)  is  the  inverse  transformation  kernel.  The  vector  x  is  related  to  x  by 
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X  =  x  +  e 


(6) 


when  e,  the  error  term,  is  0  then  the  reconstruction  is  perfect.  For  the  moment,  we  will 
require  perfect  reconstruction  and  consider  c  to  be  0  and  x  to  be  equal  to  x. 


In  matrix  notation 

y  ^  Fx 
X  =  Gy 
G  = 


for  input  vectors:,  and  coefficient  vector 3^  and 

F  =  . v-i 


we  will  use  the  notation 

G  -  . A^-1 

where  are  basis  vectors.  Therefore 

=  Gy  =  yik)gk 

k=0 


(7) 


(8) 


(9) 


X 


(10) 
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and  jr  is  the  weighted  sum  of  basis  vectors,  where  the  weights  are  the  values  of  the  transform 
coefficients. 


Let  F  be  real,  the  class  of  orthogonal  transforms  is  defined  by 

=  F^  (11) 


implying 


F^F  =  FF^  =  I 


(12) 


where  I  is  the  identity  matrix  of  order  N.  A  real  matrix  is  orthogonal  if  and  only  if  its 
columns  and  rows  form  an  orthonormal  set  where  every  vector  has  unit  length,  i.e. 


0  for  i * j 
1  for  i=j 


(13) 


where  is  the  Kronecker  delta  function  and  therefore  the  inverse  transform  is  just  the 
transpose  of  F: 


G  =  (F-y 


(14) 


Orthogonality  is  a  necessary  property  for  basis  vectors  that  are  used  to  decompose  an 
input  into  uncorrelated  components  in  an  A^-dimensional  space.  Orthonormality  is  an  even 
stronger  property;  it  leads  to  transforms  where  the  average  sum  of  the  variances  of  the 
elements  of  the  output  y  are  equal  to  the  variance  of  the  elements  of  ac.  This  implies  that 
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whatever  the  average  reconstruction  error  variance  is,  it  is  equal  the  error  variance  introduced 
in  the  quantization  of  the  transform  coefficients.  While  it  is  true  that  finite  precision  of  the 
calculations  can  introduce  errors  also,  we  will  consider  them  insignificant. 


The  two-dimensional  transform  can  be  generalized  from  (4)  and  (5) 

N-l  V-1 

yik,  0  =  S  S  ")  ")  (15) 

m=0  n=^0 


N-l  N-i 

x{m,n)  -  8(k,l  ,m,n)  (16) 

A:=0  /=0 


each  of  the  above  describe  a  element  sequence  with  transformation  filters  of  J{k,  I, m,n)  and 
g{k,l,m,n)  respectively.  For  a  2-dimensional  image  transform  we  will  generally  assume  a 
NxN  image  and  that  the  transform  filters  are  separable.  Separable  filters  allow  horizontal 
and  vertical  operations  to  be  done  separately  and  as  such  the  filters  can  be  decomposed  into 


f{k,  l,m,n)  =  fjjc,  m)  f^(l,  n) 
g(k,l,m,n)  =  g^(k,m)  g^j^Un) 


(17) 


This  permits  us  to  perform  the  transform  in  two  one-dimensional  operations.  The  two 
dimensional  transform  to  obtain  y  can  now  be  done  as 
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N-\  N-l 

y(k,l)  =  4m,«) /^(/,n) 

m=0  n=0 

N~l 

m=0 


(18) 


Defining  and  F*  as  arrays  of  the  appropriately  shifted  filters  and /*  and  requiring 
the  filters  to  be  symmetrical  allows  us  to  establish 

P  (19) 

Let  JiT  and  Y  represent  the  arrays  containing  the  elements  x(m,n)  and  y(k,l)  respectively,  and 
we  can  further  say 

Y  =  FXF^  (20) 


and  with  F'^  =  F^, 


X  =  F^YF  (21) 

Observe  that  F  is  simply  the  one-dimensional  transform  and  that  this  operation  is  only 
possible  with  separable  filters. 
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Sub-band  Transforms 

Now  with  the  understanding  of  how  the  transform  operates,  we  need  to  introduce  the 
concept  of  sub-band  transforms.  These  transforms  convert  the  original  signal  into  two  or 
more  sub-bands.  We  will  only  be  concerned  with  band  splitting  or  two  bands,  although  the 
process  can  be  applied  to  multiple  sub-bands. 

In  general,  the  band  splitting  process  involves  convolving,  or  applying  the  bandpass 
transform  to  the  input  elements,  and  then  decimating  the  resulting  output.  Decimation,  also 
known  as  sub-sampling,  samples  and  retains  only  every  ^h  element  of  the  output,  in  our  case 
we  will  use  a  two-channel  sub-band  transform  and  sub-sample  by  two.  This  means  that  we 
will  use  a  pair  of  specially  designed  filters,  a  high  bandpass  filter  G  and  low  bandpass  filter 
H.  Convolving  and  decimating  the  input  results  in  a  sub-band  output  of  NI2  elements  for  a 
N  element  input,  so  there  is  a  lowpass  output  of  NH  elements  and  a  high-pass  output  of  NI2 
elements. 

To  illustrate,  using  the  z  transform 

v-i 

x{n)  -  x(z)  =  5^  x{n)  z (22) 

n=0 


suppose  we  have  a  signal  m(0),  m(  1 ), u{2), ...,  u{2N-\)  that  we  wish  to  decimate  or  sub-sample 
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v(n)  -  u(2n)  -  -^(uin) +{-iyu(n)) 


(23) 


then  we  can  say 


N-l 

v(z)  =  Y, 

n=0 


(24) 


if 

^{u(n)  +  (-!)"«(«))  =  v(0)z‘°+0-z’'^ ' +v(l)z''^  ' +... 

+  v(A^-l)z''^  +0-v(z)‘" 

=  {(m(0)  +  (-  1)V0))z  -0  +  1(m(1)  +  (-  1)1m(1))z  '/2-'  + 

^  z 

+  ^m(2A^-1)  +(-1)2^-‘m(0))z 

=  ^M(0)fc  +  m(  1  )t  ‘^2)'  ‘  + . . .  ] 

+  1[(  - 1 )  VO)z -°  +  ( - 1 )  *  m(  1  )fe ''t  ‘  +  -  ] 

=  ^[m(z'^^)+m(-z'^^)] 

=  v(z) 

Therefore,  every  other  sample  is  used  during  sub-sampling.  The  interpolation  or  up- 
sampling  is  similar  in  approach,  using  v(0),  v(l),  v(2),...,v(A^-l)  the  sampled  signal  we  want 


to  recreate  u(0),u(l),u(2),...,u(2N-i) 
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To  reconstruct  from  the  transformed  elements,  we  must  first  upsample  or  interpolate, 
we  do  this  by  inserting  zeros  between  the  N  samples,  such  that 


m(0)  =  v(0) 

m(1)  =  0 

u{2)  =  v(l) 

m(3)  =  0  (26) 

u{2N-2)  =  v{N-\) 
u{2N-\)  =  0 


this  is  done  by 


2n-l 

m(«)  =  v(rt)z“"  =  v(0)+0(z“*)  +  v(l)z“^+(0)z"^  +  ... +v(iV-l)z^^“* 

n=0 


-  v(0)+v(l)z-2  +  ...+v(iV-l)z2^-‘ 
=  v(0)  +  v(l)z^  ^  + ... +v(A^-l)z^^  ' 

=  V(z2) 


We  can  then  convolve  the  high  and  low  pass  samples  with  reconstruction  filters  G '  dxtdH' 
respectively,  and  then  add  the  two  signals  together  to  obtain  the  original  input.  The  forward 
application  stage  of  the  sub-band  transforms  is  called  analysis  and  the  inverse  application 
is  called  synthesis  as  shown  in  Figure  1. 

The  process  of  sequentially  reprocessing  the  outputs  of  the  analysis  stage  by  another 
application  of  the  analysis  stage  is  known  as  a  uniform  cascade  system,  when  the  analysis 
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stage  is  only  reapplied  to  the  low-pass  output  of  the  analysis,  it  is  known  as  a  non-uniform 
or  pyramid  cascading  system. 


I 


Figure  1:  The  two-channel  band-pass  transform  operation. 


What  are  wavelets? 

Wavelets  are  a  family  of  functions  that  satisfy  certain  requirements.  The  name  wavelet 
should  bring  to  mind  a  function  that  waves  above  and  below  the  x-axis  and  should  therfore 
integrate  to  zero.  Another  connotation  of  wavelet  suggests  that  it  should  be  small,  or  in  other 
words  the  function  has  to  be  well  localized.  Other  constraints  on  the  functions  insure  quick 
and  easy  calculation  of  the  direct  and  inverse  wavelet  transform. 
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There  are  a  broad  variety  of  wavelets,  from  smooth  wavelets,  compactly  supported 
wavelets,  wavelets  with  simple  mathematical  expressions,  wavelets  with  simple  associated 
filters,  etc.  The  most  simple  is  the  Haar  wavelet,  and  we  discuss  it  as  an  introductory 
example.  Like  sines  and  cosines  in  Fourier  analysis,  wavelets  are  used  as  basis  functions  in 
representing  other  functions.  Once  the  wavelet  (sometimes  called  the  mother  wavelet)  ij/(x) 
is  fixed,  a  basis  can  be  made  of  translations  and  dilations  of  the  mother  wavelet 
(a,b)  e  R*x  R.  It  is  convenient  to  take  special  values  for  a  and  b  in  defining  the 
wavelet  basis:  a  -  2"^  and  £>  =  k  ■  2L  where  ^  and y  are  integers.  This  choice  of  a  and  Z? 
will  give  a  sparse  basis.  This  choice  also  connects  multiresolution  analysis  in  signal 
processing  with  the  world  of  wavelets. 

Why  not  use  the  traditional  Fourier  methods?  The  Fourier  transform  has  been  used  for 
many  years  in  signal  analysis.  It  has  proven  invaluable  in  applications  ranging  from  pattern 
recognition  to  image  processing.  There  are  certain  undesirable  limitations  that  are  inherent 
in  it  though,  and  the  wavelet  transform  has  the  same  power  and  versatility,  but  without  the 
limitations. 

The  Fourier  basis  functions  are  localized  in  frequency  but  not  in  time.  Small  frequency 
changes  in  the  Fourier  transform  will  produce  changes  everywhere  in  the  time  domain.  In 
contrast,  wavelets  are  local  in  both  frequency/scale  (via  dilations)  and  in  time  (via 


translations). 
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The  Fourier  transform  works  on  the  assumption  that  the  original  time  domain  function 
(i.e.,  signal  or  image)  is  periodic  in  nature.  As  a  result,  the  Fourier  transform  has  difficulty 
with  transient  components  of  the  image  that  are  localized  in  time.  In  other  words,  it  poorly 
localizes  sharp  transitions  in  the  image  such  as  edges.  The  Fourier  transform  of  a  signal  does 
not  convey  any  information  about  the  translation  of  the  signal  in  time.  The  only  way  to 
handle  this  shortcoming  is  to  window  the  input  data  so  that  the  sampled  data  converges  to 
zero  at  the  endpoints. 

Wavelet  transforms  provide  a  more  compact  representation  of  many  classes  of 
functions.  For  example,  functions  with  discontinuities  and  functions  with  sharp  spikes 
usually  take  substantially  fewer  wavelet  basis  functions  than  sine-cosine  basis  functions  to 
achieve  a  comparable  approximation.  This  sparse  coding  makes  wavelets  excellent  tools  in 
image  compression. 

In  recent  years,  new  families  of  orthonormal  basis  functions  have  been  developed  to 
overcome  these  shortcomings.  This  family  of  functions  are  called  ''wavelets.''  Unlike  the 
sine  and  cosine,  wave  of  the  Fourier  transform  and  discrete  cosine  transform,  they  need  not 
have  infinite  duration.  They  can  be  non-zero  for  only  a  small  range  of  the  function.  This 
short  duration  is  referred  to  as  “compact  support”  and  allows  a  time  domain  function  to  be 
translated  into  a  representation  that  is  localized  in  frequency  and  time.  This  is  the  behavior 
that  has  provided  new  advances  in  signal  processing. 
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Large  and  noisy  data  sets  can  be  easily  and  quickly  transformed  by  the  DWT.  The  data 
are  coded  by  the  wavelet  coefficients.  In  addition,  the  DWT  is  much  faster  than  the  DFT. 
It  is  well  known  that  the  computational  complexity  of  the  DFT  is  0(n  log2  n)  [25].  For  the 
discrete  wavelet  transform  the  computational  complexity  goes  down  to  0(n). 

How  do  the  wavelets  work? 

The  Haar  wavelet 

i.o 

0.5 

0.0- 

■0.5 

-1.0 

Figure  2:  Haar  Wavelet 


To  explain  how  wavelets  work,  we  start  with  an  example.  We  choose  the  simplest  and 
best  known  of  all  wavelets,  the  Haar  wavelet,  ^x).  The  Haar  wavelet  has  been  known  for 
more  than  eighty  years  and  has  been  used  in  various  mathematical  fields.  It  is  known  that 
any  continuous  function  can  be  approximated  uniformly  by  Haar  functions.  It  is  a  step 
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function  taking  values  1  and -1,  on  [0,1/2)  and  [1/2,1),  respectively.  The  graph  of  the  Haar 
wavelet  is  given  in  Figure  2. 

Dilations  and  translations  of  the  function  if/, 

iJ/jiiix)  ^  const- (28) 

where  j  is  the  dilation  function  and  k  is  the  translation  function,  define  an  orthogonal  basis 
in  L^(R)  (the  space  of  all  square  integrable  functions).  This  means  that  any  element  in  L^(R) 
may  be  represented  as  a  linear  combination  (possibly  infinite)  of  these  basis  functions. 

The  orthogonality  of  rfrji^  is  easy  to  check.  It  is  apparent  that 

j%k^x)-Hry^ix)  =  Q  (29) 

whenever^  =7 'and  k  =  k''\snot  satisfied  simultaneously.  If  j  j'  (say  j'  <  j),  then  nonzero 
values  of  the  wavelet  are  contained  in  the  set  where  the  wavelet  is  constant.  That 
makes  integral  (29)  equal  to  zero. 

Ify  =j',  but  k  9^  k',  then  at  least  one  factor  in  the  product  is  zero.  Thus  the 

functions  ifri^j  are  orthogonal.  The  constant  that  makes  this  orthogonal  basis  orthonormal  is 
2)'^.  Indeed,  from  the  definition  of  norm  in  L^: 


1 


(constfJ'f//^(2'x-k)dx  =  (constf-2'^Ji//^(t)dt  =  (const)^-2'-' 


(30) 
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Figure  3:  Dilations  and  translations  of  the  Haar  wavelet  on  [0,1] 


The  functions  ^23  shown  in  Figure  3.  The  set 

A:eZ}  defines  an  orthonormal  basis  for  L^.  Alternatively  we  will  consider 
orthonormal  bases  of  the  form  j^jg,  keZ},  where  (pg^  is  called  the  scaling 

/wncfton  associated  with  the  wavelet  basis  The  set  {0^.  keZ}  spans  the  same  subspace 

as  ^ji^,  j<jQ,  keZ}.  We  will  later  define  0;^.  For  the  Haar  wavelet  basis  the  scaling 
function  is  very  simple.  It  is  unity  on  the  interval  [0,1),  i.e. 

(p(x)  =  l(0<x<l)  (31) 


y  =  CVo’)'!’  —  2'- -i)  ^  the  data  vector  of  size  2".  The  data  vector  can  be  associated 
with  a  piecewise  constant  function /on  [0,1)  generated  by  y  as  follows. 
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2"-l 


fix)  =  E  >',-l(^2-”<x<(^+l)2-”) 

(t=0 


(32) 


The  (data)  function/is  obviously  in  the  [0, 1 )  space,  and  the  wavelet  decomposition 
of /  has  the  form: 

fix)  =  c^<Pix)  +  E  E  djM^)  (33) 

7=0  k=0 


The  sum  with  respect  to  j  is  finite  because  /is  a  step  function,  and  everything  can  be 
exactly  described  by  resolutions  up  to  the  (« -  l)-st  level.  For  each  level  the  sum  with  respect 
to  k  is  also  finite  because  the  domain  of/is  finite.  In  particular,  no  translations  of  the  scaling 
function  (pgg  are  required. 

An  obvious  disadvantage  of  the  Haar  wavelet  is  that  it  is  not  continuous,  and  therefore 
it  is  not  the  best  choice  for  representing  smooth  functions.  Therefore,  we  need  to  extend  this 
concept  to  additional  wavelet  representations. 

What  are  the  requirements  for  a  filter? 

Daubechies  discovered  that  the  wavelet  transform  can  be  implemented  with  a  specially 
designed  pair  of finite  impulse  response  (FIR)  filters  called  a  quadrature  mirror  filter  ( QMF). 
A  FIR  performs  the  dot  or  sum  of  products  between  the  filter  coefficients  and  the  discrete 
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samples  in  the  tapped  delay  line  of  the  filter.  Passing  a  set  of  discrete  samples  through  a  FIR 
filter  is  a  discrete  convolution  of  the  signal  with  the  filter’s  coefficients. 

To  picture  this  process  we  consider  a  four  coefficient  filter  designed  by  Daubechies. 
The  transform  matrix  for  this  filter  is 

Cq  c,  C2  C3  0  0  0  ...  0  0  m, 

C3  -C2  c,  -Cq  0  0  0  ...  0  0  M2 

0  0  Cq  C,  C2  C3  0  ...  0  0  M3 

0  0  C3  -C2  c,  -Cq  0  ...  0  0  M4 

(34) 

0  0  0  0  ...  0  Cq  C,  C2  C3  M^_3 

0  0  0  0  ...  0  C3  -C2  C,  -Cq  m„_2 

C2  C3  0  0  ...  0  0  0  Cq  C,  M^., 

C,  -Cq  0  0  ...  0  0  0  C3  -C2  M^ 

There  are  several  points  to  note  about  this  matrix.  The  odd  rows  are  copies  of  the  first  row 
shifted  by  two  each  time.  The  even  rows  are  not  the  same  filter,  but  do  use  the  same 
coefficients,  therefore  the  even  rows  actually  perform  a  different  convolution.  Note  that  the 
filter  wraps  around  the  matrix  boundary.  The  overall  effect  of  the  matrix  is  to  do  two 
separate  convolutions  and  to  decimate  the  results  by  two.  The  results  are  interleaved  with 
each  other.  The  filter  Cg,  c,,  Cj,  c^,  is  an  averaging  or  smoothing  filter  generally  called  H  and 
the  filter  c^,  -C2,  c ;,  -Cq,  does  not  average  because  of  the  minus  signs  and  therefore  is  the  detail 
filter  called  G. 
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There  are  some  requirements  that  must  be  placed  on  the  filter  to  make  them  perform 
properly.  First,  the  c’s  must  be  chosen  so  that  the  detail  filter  G  will  produce  a  zero  response 
to  a  sufficiently  smooth  signal.  In  the  case  of  a  four  coefficient  filter  this  means  that  we  must 
have  vanishing  moments  of  order  two.  Another  condition  is  required  to  reconstruct  the 
original  input  from  the  output  is  that  the  matrix  must  be  orthogonal.  Hence,  the 
reconstruction  matrix  for  equation  (34)  is 
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The  only  way  that  (34)  and  (35)  can  be  orthogonal  is  if  the  following  equations  are 
satisfied 


Cq  +C|  +C3  +C3 


C(jC2+C,C3  =  0 


(36) 


Additionally,  for  the  vanishing  moments  to  be  satisfied  requires  that 
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Cj  C2  +  Cj  Cg  0 
0(c3)-l(c2)+2(c,)-3(cg)  =  0 


(37) 


also  be  met.  The  equation  sets  in  (36)  and  (37)  were  first  solved  by  Daubechies  [7]  and  are  listed 
as  coefficient  set  4  in  Table  1.  As  the  number  of  coefficients  increase  by  two,  the  number 
of  orthogonality  requirements  increase  by  one. 


Figure  4:  Pyramid  output  of  the  wavelet  transform 


QMF  filters  are  unique  because  the  frequency  responses  of  the  two  filters  separate  the 
high  and  the  low  frequency  components  of  the  input  data.  A  dividing  point  is  generally 
found  between  0  and  half  the  data  sampling  (Nyquist)  frequency.  The  outputs  of  the  QMF 
filter  pair  are  decimated  by  two,  in  other  words  every  other  sample  is  discarded,  see  Figure 
1.  The  high  frequency  pass  resolves  the  details  of  the  image  and  the  low  frequency  pass 
resolves  the  averages  (smoothness)  of  the  image.  The  lowpass  portion  of  the  data  set  can 
then  be  repeatedly  fed  to  the  QMF.  Figure  4  demonstrates  how  the  input  signal  is  divided 
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into  a  smooth  (s)  and  detail  (d)  component.  The  smooth  component  is  recursively  processed 
to  create  new  smooth  and  detail  components  thereby  forming  a  pyramid. 


Figure  5:  Spectral  octave  bands. 


The  recursive  pyramiding  of  an  image  divides  the  spectrum  of  the  original  signal  into 
octaves  bands  with  successively  coarser  measurements  in  time  as  the  width  of  each  spectral 
band  decreases  in  frequency,  Figure  5.  This  has  the  effect  of  folding  most  of  the  energy 
contained  in  the  full  vector  across  the  Nyquist  frequency  into  the  lower  half  of  the  vector  as 
each  level  is  processed.  The  grey  band  reflects  the  vector  processed  for  the  next  level  of 


processing. 
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Figure  6:  Haar  wavelet  approximation  of  a  sine  wave  at  six  levels  of  resolution. 

As  an  example  of  the  process  and  the  output,  a  128-point  sine  wave  is  decomposed  by 
the  Haar  basis  set  in  Figure  6  each  level  shows  an  approximation  reconstruction  for  that 
level.  Figure  7  shows  the  actual  signal  data  set  superimposed  with  the  transform 
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coefficients.  Note  how  the  amplitudes  of  the  transform  coefficients  increase  as  the  resolution 
decreases. 


Figure  7:  128  point  sine  wave  decomposition  with  Haar  basis  set. 


Mallat  [161]  has  shown  that  the  pyramid  algorithm  can  be  applied  to  the  wavelet 
transform  by  using  the  wavelet  coefficients  as  the  filter  coefficients  of  the  QMF  filter  pairs. 
The  same  wavelet  coefficients  are  used  in  both  the  high-pass  and  lowpass  filters.  The 
lowpass  coefficients  are  associated  with  of  the  scaling  function  (p. 

(pit)  =  E  ak^i2t-k) 

k€Z 


(38) 
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The  output  of  each  lowpass  filter  is  the  approximation  components  for  that  level  of  the 
tree.  The  high-pass  filter  is  associated  with  the  of  the  wavelet  function  ifr  (note  the 
alternating  sign  change). 


keZ 


(39) 


The  output  of  the  high-pass  filter  is  the  detail  components  of  the  original  signal  at 
resolution  2^,  where  j  is  the  level  of  the  pyramid.  The  lowpass  output  of  the  previous  level 
is  used  to  generate  a  new  set  of  high-pass  and  lowpass  outputs  for  the  next  level  of  the  tree. 
Decimation  by  two  corresponds  to  the  multiresolutional  nature  (the  j  parameter)  of  the 
scaling  and  wavelet  functions. 

What  do  wavelets  look  like? 

The  following  are  visualizations  and  listings  of  the  Daubechies  filter  set.  The  filter 
visualizations  were  created  with  a  single  impulse  of  1000.0  in  the  sixth  component  a  128- 
point  data  set  in  the  frequency  domain.  This  single  impulse  was  then  returned  to  the  time 
domain,  much  in  the  same  fashion  that  pure  sine  waves  are  generated  from  a  single  impulse 
in  the  frequency  domain  and  then  returned  to  the  time  domain  with  the  Fourier  transform. 
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Figure  8:  Daubechies’  wavelet  bases  2-8. 


_ Id4  d6  d8  d10 


0.7071067811871  0.48296291 31 45 1  0.3326705529501  0.2303778133091  0.160102397974 


0.7071067811871  0.8365163037381  0.80689150931 1 1  0.7148465705531  0.603829269797 


0.2241438680421  0.459877502118  0.630880767930  0.724308528438 


-0.129409522551 1  -0.1 35011 020010  -0.027983769417  0.138428145901 


-0.085441273882  -0.187034811719  -0.242294887066 


0.035226291 8821  0.030841 381 836  -0.032244869585 


0.032883011667  0.077571493840 


-0.01 0597401 785 1  -0.006241 49021 3 


-0.012580751999 


0.003335725285 


Table  1:  Daubechies’  wavelet  coefficients  2-10. 


16  32  48  64  80  96  112  128 


0  16  32  48  64  80  96  112  128 


Figure  9:  Daubechies’  wavelet  bases  10,  12,  16  and  20. 


d12 


0.111540743350 


0.494623890398 


0.751133908021 


0.315250351709 


-0.226264693965 


-0.129766867567 


0.097501605587 


0.027522865530 


-0.031582039318 


0.000553842201 


0.004777257511 


-0.001077301085 


0,077852054085 


0.396539319482 


0.729132090846 


0.469782287405 


-0.143906003929 


-0.224036184994 


0.071309219267 


0.080612609151 


-0.038029936935 


-0.016574541631 


0.012550998556 


0.000429577973 


■0,001801640704 


0.000353713800 


d16 


0.054415842243 


0.312871590914 


0.675630736297 


0.585354683654 


-0.015829105256 


-0.284015542962 


0.000472484574 


0.128747426620 


-0.017369301002 


-0.044088253931 


0.013981027917 


0.008746094047 


-0.004870352993 


-0.000391740373 


0.000675449406 


-0.000117476784 


d18 


0.038077947364 


0.243834674613 


0.604823123690 


0.657288078051 


0.133197385825 


-0.293273783279 


-0.096840783223 


0.148540749338 


0.030725681479 


-0.067632829061 


0.000250947115 


0.022361662124 


-0.004723204758 


-0.004281503682 


0.001847646883 


0.000230385764 


d20 


0.026670057901 


0.188176800078 


0.527201188932 


0.688459039454 


0.281172343661 


-0,249846424327 


-0.195946274377 


0,127369340336 


0.093057364604 


-0.071394147166 


-0.029457536822 


0.033212674059 


0.003606553567 


-0.010733175483 


0.001395351747 


0.001992405295 


0.000251 9631 891-0.000685856695 


-0.000116466855 


0.000093588670 


-0.000013264203 


Table  2:  Daubechies’  wavelet  coefficients  12-20. 
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Quantization 

Quantization  is  the  process  of  converting  a  given  symbol  value  into  a  value  taken  from 
a  set  of  finite  possibilities.  Up  to  this  point  all  of  the  methods  discussed  are  lossless  except 
for  finite  precision  arithmetic  errors.  It  is  normally  the  step  in  a  compression  scheme  that 
introduces  the  noise  (loss)  into  the  system.  Consider  as  a  form  of  quantization,  the 
conversion  from  real  numbers  to  integers,  the  integers  are  finite  for  a  given  range  and  the  loss 
of  the  fractional  part  of  the  real  number  is  the  noise  that  is  introduced.  Quantization  allows 
the  number  of  symbols  to  be  reduced  and  correspondingly  improves  encoding  efficiency. 

After  using  the  DWT  transform  on  an  image  the  resulting  array  will  be  in  floating  point 
representation.  We  have  quadrupled  the  size  of  the  image  since  the  floating  point 
representation  requires  in  general  four  times  the  storage  (four  bytes  output  versus  one  byte 
input).  Now,  even  to  represent  the  image  in  the  same  space  as  the  original  will  require 
conversion  to  a  finite  precision  representation.  It  is  at  this  point  that  we  begin  to  introduce 
noise  into  the  image  and  a  corresponding  loss  of  fidelity. 

First,  let  us  assume  that  the  output  range  of  our  transform  can  be  truncated  or  rounded 
to  fit  into  an  integer  representation.  Using  our  model  from  before,  the  quantization  step  is 
inserted  between  transform  steps  as  in  Figure  10. 
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■  I 


Figure  10:  Transform  model  with  quantization  step. 

One  of  the  key  features  of  the  wavelet  transform  is  that  the  coefficients  can  be  severely 
truncated  and  still  reconstruct  an  adequate  representation  of  the  input.  The  next  step  in  this 
process  then  is  to  truncate  or  threshold  the  data  so  that  a  minimal  set  is  retained. 

Maximizing  Retained  Data 

It  is  in  this  portion  of  the  scheme  that  we  wish  to  address  the  problem  of  determining 
which  coefficients  we  wish  to  keep  to  represent  the  image  for  an  expected  byte  budget. 
Recall  for  a  moment  that  the  transform  removes  the  smoothness  of  the  image  at  each  level 
of  the  transform  and  its  inability  to  absorb  all  of  the  changes  in  the  image  are  the  coefficients 
that  are  left  behind  in  the  high  pass  portion  of  the  sample  set.  If  the  details  are  small  enough 
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then  they  can  be  omitted  without  significantly  affecting  the  overall  quality  of  the 
reconstructed  data. 

This  omission  of  small  details  is  called  thresholding.  In  general,  one  of  two  methods 
are  used,  hard  thresholding  is  when  any  coefficient  value  below  an  established  threshold  is 
reduced  to  zero.  Soft  thresholding  is  a  shrinking  of  values  towards  zero,  any  coefficient 
below  a  certain  value  is  reduced  to  zero  and  all  other  values  are  shifted  towards  zero  by  the 
threshold  amount.  The  threshold  can  be  set  by  a  fixed  value  or  by  a  percentile  value  based 
on  the  total  range  of  the  coefficients.  For  example,  if  we  wish  to  reduce  the  number  of 
coefficients  by  half  then  we  set  the  threshold  at  the  median  value.  It  is  thresholding  along 
with  the  finite  representation  from  quantization  that  will  allow  us  to  further  compress  our 
image. 

From  our  discussion  of  the  transform,  each  level  of  depth  of  the  transform  is  at  a 
coarser  resolution  by  a  factor  of  two.  The  detail  coefficients  in  the  highband  pass  of  a  given 
level  are  only  correlated  to  the  highband  pass  of  the  next  level  of  the  transform  in  the  sense 
that  if  a  level  has  an  area  of  significant  detail  coefficient  activity  then  it  indicates  that  there 
may  be  more  significant  detail  coefficients  at  a  finer  resolution  level.  Accordingly,  if  the 
coefficients  are  small  then  the  detail  probably  contributes  little  to  a  relatively  smooth  area 
of  the  image.  What  this  means  is  simply,  if  there  are  a  lot  of  edges  or  variations  in  an  area 
at  a  given  resolution  then  it  is  a  good  indicator  of  additional  detail  at  a  finer  scale. 
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Conversely,  smoothly  transitioning  areas  will  not  gain  significant  detail  at  greater  levels  of 
resolution.  We  will  refer  to  this  concept  as  self -similarity. 

We  will  exploit  the  self-similarity  between  the  various  levels  of  the  transformation. 
Although  the  transform  removes  most  of  the  correlation  of  the  image,  it  does  not  remove  all 
of  it  and  the  correlation  left  is  not  necessarily  predictable.  It  is  this  self-similar  relationship 
between  levels  we  will  use  to  model  our  storage  of  the  image.  The  concept  is  straight 
forward,  use  a  binary  tree  to  represent  each  level  of  the  transformation  with  the  parent 
pointing  to  the  children  that  are  the  coefficients  at  the  next  level  of  higher  resolution.  This 
is  similar  to  research  done  by  J.  M.  Shapiro  [30]  except  that  they  used  a  quadtree  model  with 
only  a  single  decision  level.  If  the  value  of  the  coefficient  is  significant,  then  keep  it.  If  the 
value  of  the  coefficient  is  not  significant,  nor  are  any  of  the  values  of  the  descendants,  then 
discard  it  and  prune  the  branch.  Thus,  Figure  11  shows  our  modified  compression  scheme 
model. 


Figure  11;  Model  with  tree  storage  step. 
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Initially,  the  problem  is  how  to  build  the  decision  process  into  the  tree  structure. 
Special  symbols  can  be  allocated  from  the  available  alphabet  to  direct  the  construction  and 
reconstruction  of  the  binary  tree.  Alternatively,  a  bit  prefix  can  be  used  from  the  bits 
allocated  for  the  numerical  representation  of  the  coefficient. 

We  will  evaluate  these  models  on  a  simple  example  using  a  four-level  transform  of  3 1 
elements.  Placing  the  root  at  the  coarsest  level  (deepest  level  of  the  transformation),  we  can 
imagine  that  the  same  relative  positions  in  the  finer  levels  (earlier  levels  of  the 
transformation),  create  a  binary  tree  structure.  Thus,  34  is  the  lowest  level  of  the 
transformation  and  the  root  of  the  points  above  it. 
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Figure  12:  Four  level  transformation  map. 


We  will  set  the  thresholding  level  at  10  and  ignore  elements  of  a  lesser  value. 
Normally,  this  would  be  an  absolute  value  threshold  since  the  coefficients  will  range  above 
and  below  zero,  but  for  simplicity  we  will  use  only  positive  values. 
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Figure  13:  Map  after  setting  threshold  at  10. 


If  we  use  stop  codes  to  determine  the  tree  structure  then  for  every  non-terminal  branch 
node,  we  must  keep  one  extra  node.  For  every  terminal  node  (except  for  the  last  level  which 
we  can  predetermine)  we  must  keep  two  extra  nodes.  The  worst  case  for  this  model  is  a  full 
binary  tree  that  does  not  reach  the  last  level.  There  will  be  n  -i- 1  stop  codes  for  n  nodes  kept. 
This  clearly  makes  the  overhead  greater  than  50%  in  the  worst  case.  The  best  case  would  be 
a  full  binary  tree  where  every  node  is  kept  and  since  no  stop  codes  would  be  necessary  the 
overhead  would  be  zero.  Since  the  purpose  of  the  tree  encoding  is  to  reduce  the  number  of 
nodes  stored  then  this  is  non-optimal. 


34 

22 

25 

15 

18 

22 

10 

12 

■ 

4 

16 

12 

18 

■ 

■ 

□D 

11  ■ 

10  12 

■  11 

13  14 

Figure  14:  Stop  codes  for  four  level  transformation  map. 


The  choice  then  is  how  many  bits  to  allocate  to  the  decision  process.  If  only  one  bit 
is  allocated,  then  it  becomes  either  a  “go/no”  decision  on  whether  to  pursue  subsequent 
descendent  levels.  If  the  decision  is  to  go  because  one  or  more  of  the  child  nodes  are 
significant,  then  we  must  keep  both  child  nodes  even  if  only  one  is  significant.  If  we  allow 
two  bits  for  the  decision  we  will  half  the  remaining  precision  left  in  the  byte  storage,  64 
possibilities  instead  of  128,  but  we  will  improve  the  overhead  of  the  decision  process  by 
reducing  the  number  of  subsequent  bytes  that  must  be  kept,  see  Table  3. 
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Figure  15:  Map  with  1-bit  binary  tree  decision  prefix. 
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Figure  16:  Map  with  2-bit  binary  tree  decision  prefix. 
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In  the  one  bit  model,  if  the  decision  is  made  to  stop  at  a  given  point  then  we  have  seven 
remaining  bits  in  the  element  for  the  value.  If  we  must  continue  to  the  next  level  because 
either  the  left,  right,  or  both  subbranches  have  significant  coefficients  then  both  subbranch 
elements  must  be  kept  at  an  additional  cost  of  two  bytes.  If  the  two  bit  model  is  used  then 
zero,  one,  or  two  additional  elements  are  kept  according  to  the  decision.  The  average 
overhead  for  the  one  bit  model  is  20  bits  while  the  average  overhead  for  the  two  bit  model 
is  only  16  bits.  Accordingly,  the  scheme  is  25%  more  efficient  if  the  resolution  of  six  bits 
instead  of  seven  can  be  used. 


Reconstruction  is  very  simple,  read  the  prefix,  if  there  is  a  1  in  the  first  bit  then 
continue  recursively  building  the  left  subtree  until  a  leftmost  0  is  encountered.  As  0  bits  are 


encountered,  go  right  (the  second  bit  is  a  1)  or  return  a  level  and  resume  the  process.  The 
depth  of  the  reconstruction  process  is  self-limiting  and  the  process  can  be  terminated  by 
exhausting  the  input. 

So  [Oil..],  would  represent  pruning  the  left  branch  and  continuing  the  right  branch. 


...  [11134]  [11122]  [10115]  [00112]  [1 1118  ][10I4  ]  [00111]  [11116]  [00110] 
[01112]  [10125]  [1 1122  ]  [00112]  [00111]  [1 1118]  [00113]  [00114]  [00110  ] ... 
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Figure  17:  Reconstructed  map  from  coded  stream. 


Reconstruction  is  very  simple,  read  the  prefix,  if  there  is  a  1  in  the  first  bit  then 
continue  recursively  building  the  left  subtree  until  a  leftmost  0  is  encountered.  As  0  bits  are 
encountered,  go  right  (the  second  bit  is  a  1 )  or  return  a  level  and  resume  the  process.  The 
depth  of  the  reconstruction  process  is  self-limiting  and  the  process  can  be  terminated  by 
exhausting  the  input. 


Chapter  4 


Color  Images 


Up  to  this  point  we  have  been  discussing  concepts  related  to  8-bit  monochrome 
greyscale  images.  The  coding  of  24-bit  color  images  is  more  complex  issue.  At  first  glance 
it  would  appear  that  we  will  need  three  times  the  efficiency  to  code  the  three  input  channels 
(R,G,B)-  In  fact,  there  is  a  great  deal  of  redundancy  in  the  three  color  planes  and  we  can 
reduce  the  redundancy  by  performing  special  transforms  on  them.  A  color  signal  can  be 
transformed  into  a  luminance  signal  and  two  color  sub-bands  as  is  done  in  television 
broadcasting.  The  bandwidth  of  the  two  color  sub-bands  can  be  much  smaller  than  the 
luminance  bandwidth  to  reproduce  the  spatial  detail  accurately,  perhaps  as  small  as  20%[3]. 
We  can  use  the  reduction  in  the  bandwidth  to  our  advantage,  concentrating  our  best 
reconstruction  efforts  on  the  luminance  signal  while  devoting  a  much  smaller  byte  budget 
to  the  color  sub-bands. 

Two  transforms  are  in  common  use  today,  those  of  the  NTSC  (U.S.  television 
standard)  and  the  PAL  (European  standard)  transforms.  The  NTSC  YIQ  transform  is  given 
by 
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Y 

I 

= 

Q 

0.299  0.587  0.114 

0.596  -0.274  -0.322 
0.211  -0.523  0.312 


(40) 


where  Y  is  the  monochrome  luminance  channel  and  the  I  and  Q  channels  carry  the  color 
opposition  signals.  The  inverse  transform  is 


1 

0.956 

0.621 

1 

-0.273 

-0.647 

1 

-1.104 

1.701 

(41) 


The  PAL  transform  is  similar  in  that  the  Y  is  derived  in  the  same  manner,  but  the  color 
difference  signals  are  directly  derived  as  follows 

f/^  =0.493(5  -  T) 

y,  =  0.877(5 -y)  .42) 


the  scale  factors  are  introduced  to  limit  the  maximum  instantaneous  signal  amplitude. 


An  alternative  transform  corresponding  to  the  CIE  uniform  chromaticity  scale  may  also 
be  applied,  where  the  V  channel  corresponds  to  the  Y  channel  of  the  YIQ  transform. 


U 

0.405  0.116  0.133 

R 

V 

= 

0.299  0.587  0.114 

G 

W 

0.145  0.827  0.627 

B 

With  a  further  conversion  of 


50 


u  = 


U 

u+v+w 


and 


U  +  V+W 

to  give  a  y,M,v  coordinate  system. 


(44) 


(45) 


Another  transform  based  on  a  KLT  coordinate  rotation  to  generate  uncorrelated  color 
components  was  developed  by  W.  K.  Pratt  [24]  and  is  given  by  the  orthogonal  transform 


0.575  0.615  0.540 

R 

- 

0.608  0.120  -0.785 

G 

^3 

0.548  -0.779  0.305 

B 

(46) 


where  the  Kj  channel  carries  the  luminance  and  the  fL  channels  carry  the  color 

opposition  components.  Its  inverse  is  given  by 


R 

0.575  0.608  0.547 

G 

= 

0.615  0.120  -0.779 

K, 

B 

0.539  -0.785  0.305 

K, 

(47) 


From  Pratt’s  work  it  was  found  that  for  an  example  image  the  signal  energy 


redistribution  is  that  shown  in  Table  4 
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System 

Component 

1  (%)  2  (%) 

3  (%) 

RGB 

33.2 

36.2 

30.6 

YIQ 

93.0 

5.3 

1.7 

K,  K3 

94.0 

5.1 

0.9 

Table  4:  Color  signal  energy  redistribution 


As  is  readily  apparent  from  the  table,  the  dynamic  ranges  of  the  color  opposition  signals  of 
either  transform  is  much  smaller  than  the  luminance  signal. 

An  example  of  the  RGB  color  planes  and  the  Kj,  Kj,  Kj  component  planes  from  a  24- 
bit  image  of  fruit  are  in  Figure  18  and  Figure  19.  In  the  RGB  images  it  is  easy  to  see  how 
much  redundant  information  is  contained  in  all  three  subimages.  In  the  Ki,  I^,  and 
component  images  the  K,  plane  contains  the  bulk  of  the  image  information,  while  the  K  jand 
Kj  channels  contain  minimal  information. 

It  was  found  that  the  Kj  and  Kj  color  opposition  channels  can  be  compressed  as  much 
as  0.063  bits/pixel  and  still  reconstruct  a  visually  reasonable  color  map  for  the  image.  It  is 
this  ability  to  restore  an  adequate  color  map  that  allows  visual  degradation  in  that  is 
unacceptable  in  a  monochrome  image  to  be  surpassed  by  a  color  reconstruction. 
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K1  plane 


Figure  19:  K,  K2  K3  component  separation  of  fruit. 
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Color  transformation  can  be  improved  by  looking  at  human  visual  systems  modeling, 
a  logarithmic  function/bandpass  of  the  form 


T, 

0.299  0.587  0.114 

R 

T2 

= 

0.607  0.174  0.201 

G 

T, 

0  0.066  1.117 

B 

and 


(48) 


21.5  0  0 

log  T, 

= 

-41.0  41.0  0 

log  T, 

.^3. 

-6.27  0  6.27 

log  T, 

(49) 


is  more  effective  at  converting  the  color  system  into  one  that  responds  to  detectable  changes 
observed  by  the  eye. 

The  T I  channel  carries  the  luminance  signal  and  the  and  Tj  channels  carry  the  color 
opposition  signals  and  are  invariant  to  linear  scaling  of  the  R,  G,  and  B  inputs.  The  Gj  and 
G}  correspond  to  the  red-green  and  blue-yellow  channels.  A  key  feature  of  this  model  is  that 
it  is  very  effective  in  compacting  the  total  R,  G,  and  B  color  energy  selectively  into  the  G, 
channel  while  leaving  the  other  two  channels  significantly  correlated. 


Chapter  5 


A  Multi-Level  Compression  Model 

Implementation 

To  test  and  develop  the  color  image  compression  model  we  implemented  the 
components  in  separate  modules  much  as  we  have  described  the  development  of  this  paper, 
see  Figure  20.  It  consists  of  four  separate  stages:  The  color  transformation  stage,  dividing 
the  image  into  three  separate  sub-planes.  The  wavelet  transform  stage,  transforming  each 
of  the  sub-planes  with  a  wavelet  transform.  Next  is  the  quantization  step,  converting  the 
transformed  planes  into  a  finite  representation.  The  last  step  is  to  compress  the 
representation  by  reducing  the  number  of  elements  stored.  Reconstruction  is  accomplished 
by  reversing  the  sequence  and  using  the  inverse  operators. 

We  elected  to  use  Pratt’s  color  transform  in  our  implementation  because  of  its 
simplicity,  good  energy  redistribution,  and  orthogonality.  The  original  24-bit  image  is 
converted  into  three  planes  of  floating  point  values.  The  transformation  is  lossless  and 
perfect  reconstruction  can  be  obtained  from  these  three  planes. 
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Figure  20:  Multi-level  color  compression  model. 


More  than  90%  of  the  total  amount  of  information  is  contained  in  the  first  or  Kj  plane, 
dictating  the  luminance  and  therefore  the  edges  in  the  image.  Hence,  it  is  from  this  plane  that 
we  must  retain  the  greatest  information.  If  a  general  mapping  of  the  color  opposition  planes 
can  be  obtained,  we  will  be  able  to  approximate  the  locations  and  variations  of  the  color. 
This  allows  us  to  sacrifice  detail  information  contained  in  these  planes.  Approximately  10%- 
20%  of  the  data  required  for  reconstruction  of  the  Kj  plane  is  required  for  the  ^  and 
planes.  Care  must  be  taken  to  ensure  that  the  same  level  of  approximation  is  used  on  Kj  and 
K3  planes  since  the  balance  of  color  opposition  must  be  maintained  to  reconstruct  a  color 
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scheme  that  is  in  the  correct  range.  Unbalancing  of  these  planes  will  cause  undesirable  color 
shifts. 


Figure  21:  Snaking  the  image. 

A  primary  consideration  for  our  model  was  flexibility  and  the  ability  to  handle  various 
image  sizes  along  with  reducing  the  implemention  complexity.  Since  two  dimensional  QMF 
models  require  images  to  be  in  N'  formats  where  N  is  some  power  of  two,  we  chose  to  use 
a  separable  linear  transform  method  and  modified  it  to  handle  various  image  formats.  The 
image  was  treated  as  a  one  dimensional  array  rather  than  a  two  dimensional  array.  The 
advantages  are  that  any  size  image  can  be  processed  although  it  can  only  be  processed  to  the 
greatest  power  of  two  that  will  span  the  image  or  it  must  be  appended  with  an  array  of  zeros 
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to  reach  the  next  highest  power  of  two.  Usually,  the  passes  of  the  filter  are  made  across  the 
image  boundary,  the  end  of  a  scan  line  to  the  next  scan  line,  this  discontinuity  of  the  image 
creates  spikes  in  the  transform  domain.  These  spikes  can  be  misidentified  as  significant 
coefficients.  To  eliminate  this  edge  effect,  an  option  to  snake  the  image  was  implemented. 
Snaking  the  image,  as  is  implied  in  Figure  21,  takes  every  other  scan  line  and  reverses  the 
order.  The  image  edges  are  now  next  to  each  other  creating  continuity  in  the  image  and  the 
only  point  that  the  boundary  condition  exists  is  in  the  wrap  around  from  the  end  of  the  image 
to  the  first  element  of  the  image.  This  effect  can  be  seen  on  the  image  “Lenna”,  Figure  22 
and  Figure  23.  The  wrap  around  is  necessary  because  as  the  filter  approaches  the  end  of  the 
image  it  must  have  sufficient  elements  to  calculate  the  coefficients.  The  over  run  is 
accomplished  by  wrapping  around  to  the  beginning  of  the  image  instead  of  adding  unwanted 
buffering  elements.  We  implemented  a  range  of  filters  with  which  to  experiment.  These  are 
the  Daubechies  filters  of  4  -  20  coefficients. 
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Quantizing  is  the  weakest  portion  of  our  model  since  it  is  only  a  linear  bin  distribution 
based  on  the  construction  of  a  set  of  quantizing  bins  from  a  fixed  number  of  bins  for  the 
given  sub-band,  see  Table  5.  Since  each  deeper  level  of  the  pyramid  contains  a  greater 
portion  of  the  energy  of  the  image  the  number  of  bins  is  increased  for  each  level  of  the 
transform.  This  captures  about  93%  of  the  values  within  the  bins  specified  for  the  level. 

This  gives  reasonable  flexibility  to  the  coefficient  distribution  and  reconstruction.  The 
total  cost  for  the  extended  coding  is  two  and  three  bytes  respectively.  Example  runs  can  be 
seen  in  Table  6  and  Table  7. 


level  16,  sect  1  int  quantizing,  100  bins 
mean  =  17008.0,  var  =  725523185,89,  stddev  =  26935.54 
high  =  6361 4.39,  low  =  -31 7.83,  size  =  4,  bins  =  1 00 

binsize  =  808.07,  zeros  =  0,  longest  zero  run  =  0 
outsiderthresh  =  80806.61,  outsiders  =  0 
For  level  16,  sect  2  quantizing,  100  bins 
mean  =  764.5,  var  =  575880.96,  stddev  =  758.87 
high  =  1744.64,  low  =  -375.55,  size  =  4,  bins  =  100 

binsize  =  22.77,  zeros  =  0,  longest  zero  run  =  0 
outsiderthresh  =  2276.60,  outsiders  =  0 
For  level  15,  sect  2  quantizing,  93  bins 
mean  =  17.1,  var  =  510768.84,  stddev  =  714,68 
high  =  1032.83,  low  =  -1 109.04,  size  =  8,  bins  =  93 

binsize  =  23.05,  zeros  =  0,  longest  zero  run  =  0 
outsiderthresh  =  2144.04,  outsiders  =  0 
For  level  14,  sect  2  quantizing,  87  bins 
mean  =  3.6,  var  =  93784.47,  stddev  =  306.24 
high  =  545.63,  low  =  -581.80,  size  =  16,  bins  =  87 

binsize  =  10.56,  zeros  =  0,  longest  zero  run  =  0 
outsiderthresh  =  918.73,  outsiders  =  0 

For  level  6,  sect  2  quantizing,  37  bins 
mean  =  1 .1 ,  var  =  32306.60,  stddev  =  1 79.74 
high  =  662.1 8,  low  =  -659.70,  size  =  4096,  bins  =  37 

binsize  =  14.57,  zeros  =  293,  longest  zero  run  =  3 
outsiderthresh  =  539.22,  outsiders  =  19 
For  level  5,  sect  2  quantizing,  31  bins 
mean=  0.4,  var  =  9105.61,  stddev  =  95.42 
high  =  390.33,  low  =  -372.62,  size  =  8192,  bins  =  31 

binsize  =  9.23,  zeros  =  941 ,  longest  zero  run  =  6 
outsiderthresh  =  286.27,  outsiders  =  54 
For  level  4,  sect  2  quantizing,  25  bins 
mean=  0.0,  var  =  2898.01,  stddev  =  53.83 
high  =  317.94,  low  =  -290.04,  size  =  16384,  bins  =  25 

binsize  =  6.46,  zeros  =  3101,  longest  zero  run  =  15 
outsiderthresh  =  161.50,  outsiders  =  238 
For  level  3,  sect  2  quantizing,  1 8  bins 
mean  =  0.3,  var  =  731.22,  stddev  =  27.04 
high  =  181.01,  low  =  -170.25,  size  =  32768,  bins  =  18 
binsize  =  4.51 ,  zeros  =  11975,  longest  zero  run  =  32 
outsiderthresh  =  81.12,  outsiders  =  699 
For  level  2,  sect  2  quantizing,  12  bins 
mean  =  0.0,  var  =  169.55,  stddev  =  13.02 
.  high  =  1 32. 1 4,  low  =  - 1 70. 1 9,  size  =  65536,  bins  =  1 2 
binsize  =  3.26,  zeros  =  36851 ,  longest  zero  run  =  50 
outsiderthresh  =  39.06,  outsiders  =  1698 
For  level  1 ,  sect  2  quantizing,  6  bins 
mean  =  0.0,  var  =  26.89,  stddev  =  5.19 
high  =  60.81 ,  low  =  -64.59,  size  =  1 31072,  bins  =  6 
binsize  =  2.59,  zeros  =  86044,  longest  zero  run  =  72 
outsiderthresh  =  15.56,  outsiders  =  2823 


Table  5:  Example  run  from  quantizer. 
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The  purpose  of  selecting  these  simple  model  parameters  is  to  explore  the  model’s 
binary  tree  storage  capabilities.  Exception  encoding  is  employed  for  values  outside  this 
range  by  encoding  a  special  symbol  and  adding  an  additional  one  or  two  b3^es  to  the  storage. 
In  this  manner  we  could  relax  the  bin  partitioning  requirements  and  allow  a  small  percentage 
of  the  quantized  values  to  exceed  the  normal  range.  The  first  special  code  allows  one  byte 
storage  of  values  up  to  ±127.  The  second  special  code  provides  a  range  of  ±32768. 


levels  =  15,  direction  1,  lookahead  =  10,  thresh  1  =  15,  thresh2  =  0 
saving  header  of  size  34 
setting  up  tree  pointers 
maxlevel  is  1 8,  levels  requested  are  1 5 
loc  1 8  offset  0  contains  2  floats 
loc  1 7  offset  2  contains  2  floats 
loc  1 6  offset  4  contains  4  fioats 
loc  1 5  offset  8  contains  8  floats 
loc  1 4  offset  1 6  contains  1 6  floats 
loc  1 3  offset  32  contains  32  floats 
loc  1 2  offset  64  contains  64  floats 
loc  1 1  offset  1 28  contains  1 28  floats 
loc  1 0  offset  256  contains  256  floats 
loc  9  offset  51 2  contains  512  floats 
loc  8  offset  1 024  contains  1 024  floats 
loc  7  offset  2048  contains  2048  floats 
loc  6  offset  4096  contains  4096  floats 
loc  5  offset  81 92  contains  8192  floats 
loc  4  offset  16384  contains  16384  floats 
loc  3  offset  32768  contains  32768  floats 
loc  2  offset  65536  contains  65536  floats 
loc  1  offset  1 31 072  contains  1 31 072  floats 
8  tokens  copied,  8  tokens  to  process 
level  15,  left  at  0,  i=-10,  o=54 
level  14,  left  at  0,  i=51,  o=51 
level  13,  left  at  0,  i=-4,  o=60 
level  12,  left  at  0,  l=-16,  o=48 
level  11,  left  at  0,  i=-6,  o=58 
level  10,  left  at  0,  i=-23,  o=41 
level  9,  left  at  0,  i=0,  o=0 
level  8,  left  at  0,  i=7,  o=7 

level  7,  right  at  0,  i=-7,  o=-71  -  pruned  left,  going  right 
level  6,  right  at  1,  i=0,  o=-128  -  pruned  left,  going  right 
level  5,  stop  at  3,  i=-15,  o=-15  -  both  pruned,  going  back 
level  8,  right  at  0 

level  7,  stop  at  1 ,  i=-4,  o=-4  -  both  pruned,  going  back 

level  9,  right  at  0 

level  8,  left  at  1,  i=-3,  o=61 

level  7,  left  at  2,  i=2,  o=2 

level  6,  left  at  4,  i=0,  o=64  -  pruned  right,  going  left 
level  5,  right  at  8,  i=-3,  o=-67  -  pruned  left,  going  right 
level  4,  stop  at  17,  i=16,  o=-48  -  both  pruned,  going  back 
level  6,  right  at  4  skipped 

output  file  write  of  25871 

Ienna.tki.q2  =  262144  tokens,  converted  to  25871  char 
0  long  codes  output 

output  stop  =  8608,  left  =  4041 ,  right  =  391 3,  spec  =  385 
total  code  tokens  1 6947 
special  token  overhead  1 4.8% 


Table  6:  Binary  tree  encoding  output. 


levels  =  15,  direction  -1, saving  header  of  size  34 
reading  in  25984  data  tokens 
rebuilding  tree  pointers 
maxlevel  is  18,  levels  to  restore  are  15 
loc  18  offset  0  contains  2  floats 
loc  17  offset  2  contains  2  floats 
loc  1 6  offset  4  contains  4  floats 
loc  1 5  offset  8  contains  8  floats 
loc  14  offset  16  contains  16  floats 
loc  1 3  offset  32  contains  32  floats 
loc  1 2  offset  64  contains  64  floats 
loc  1 1  offset  1 28  contains  1 28  floats 
loc  10  offset  256  contains  256  floats 
loc  9  offset  512  contains  51 2  floats 
loc  8  offset  1 024  contains  1 024  floats 
loc  7  offset  2048  contains  2048  floats 
loc  6  offset  4096  contains  4096  floats 
loc  5  offset  81 92  contains  81 92  floats 
loc  4  offset  1 6384  contains  1 6384  floats 
loc  3  offset  32768  contains  32768  floats 
loc  2  offset  65536  contains  65536  floats 
loc  1  offset  1 31 072  contains  1 31 072  floats 
8  tokens  copied 
level  15,  left  at  0,  i=54,  o=-10 
level  14,  left  at  0,  i=51,o=51 
level  13,  left  at  0,  i=60,  o=-4 
level  12,  left  at  0,  i=48,  o=-16 
level  1 1 ,  left  at  0,  i=58,  o=-6 
level  10,  left  at  0,  i=41,  o=-23 
level  9,  left  at  0,  i=0,  o=0 
level  8,  left  at  0,  i=7,  o=7 

level  7,  right  at  0,  i=-71,  o=-7  -  pruned  left,  going  right 
level  6,  right  at  1,  i=-128,  o=0  -  pruned  left,  going  right 
level  5,  stop  at  3,  i=-15,  o=-15  -  both  pruned,  going  back 
level  8,  right  at  0 

level  7,  stop  at  1,  i=-4,  o=-4  -  both  pruned,  going  back 

level  9,  right  at  0 

level  8,  left  at  1 ,  i=61 ,  o=-3 

level  7,  left  at  2,  i=2,  o=2 

level  6,  left  at  4,  i=64,  o=0  -  prune  right,  going  left 

level  5,  right  at  8,  i=-67,  o=-3  -  pruned  left,  going  right 

level  4,  stop  at  17,  i=-48,  o=16  -  both  pruned,  going  back 

level  6,  right  at  4  skipped 

level  15  ,size  8,  index  8  completed 

Ienna.tki.q2.zt2  =  25984  tokens,  converted  to  25486  floats 

output  file  write  of  262144 


Table  7:  Binary  tree  reconstruction  output. 
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A  visual  example  of  the  binary  tracing  of  the  coefficients  is  in  Figure  24.  When 
looking  at  this  example  keep  in  mind  that  each  level  is  actually  one  half  the  size  of  the 
previous  level.  In  the  Lenna  example.  Figure  25,  the  image  width  spreads  these  blocks 
across  the  original  width  of  the  image.  Tracings  are  used  to  give  an  idea  of  the  relationship 
of  the  blocks  to  each  other. 


Figure  24:  Binary  relationship  between  blocks. 
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Figure  25:  Example  coefficient  binary  tree  relationship  tracings. 


An  example  of  the  restored  image  of  lenna  after  8.8: 1  compression  is  in  Figure  26. 
At  this  level  of  compression  the  losses  are  very  minor.  Note  the  loss  in  the  area  around  the 
top  of  the  hat  along  with  some  edge  degradation  around  the  frame  of  the  mirror  and  the 
feathers  of  the  hat  as  compared  to  the  image  in  Figure  22. 
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Figure  26:  Restored  image  of  Lenna  after  8.8;  1  compression. 


We  evaluated  the  performance  of  the  color  transforms  and  the  K,,  K2,  K3  transform 
provided  the  best  performance  among  it  and  the  CIE,  YIQ,  and  PAL  transforms.  Because 
of  the  crudeness  of  the  binary  tree  model  it  was  necessary  to  switch  out  the  compression 
engine  to  effectively  evaluate  the  color  transformations.  By  changing  the  compression 
engine  to  a  two-dimensional  run-length  model  compression  ratios  were  significantly 
enhanced.  We  found  the  K,,  K2,  model  to  be  the  least  sensitive  to  color  component 
degradation  as  long  as  the  K2  and  K3  channels  are  compressed  with  the  same  parameters. 
The  K3  channel  generally  required  only  half  the  number  of  coefficients  as  K2channel.  Since 
the  K|  channel  carries  more  than  90%  of  the  information,  the  other  two  channels  can  be 
compressed  to  less  than  half  of  the  K,  channel. 
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Performance 

The  binary  tree  storage  proved  to  be  a  viable  concept,  compression  ratios  of  10: 1  for 
8-bit  images  were  achieved  even  though  the  model  only  used  a  horizontal  pass  for  the 
pyramid  construction.  With  the  addition  of  horizontal  and  vertical  passes,  energy  movement 
should  be  significantly  greater  and  the  coefficient  reduction  in  the  highband  pass  should  be 
closer  to  zero.  The  most  difficult  process  is  the  discrimination  of  significant  coefficients 
between  levels  of  the  pyramid.  Since  the  energy  is  packed  into  fewer  coefficients  at  each 
level,  the  coefficients  in  the  lower  levels  of  the  pyramid  contribute  more  to  the  overall  energy 
of  the  image  than  those  of  the  higher  levels.  Care  must  be  used  when  eliminating 
coefficients  at  the  lower  levels  because  the  removal  can  cause  significant  dropouts  in  the 
affected  portions  of  the  image.  Conversely,  the  coefficients  at  the  higher  levels  establish  the 
details  of  the  image  and  must  be  carefully  retained  to  prevent  washout  of  the  image. 

Although  the  quantizer  used  is  very  crude,  it  was  found  that  by  first  determining 
significance  across  the  levels  of  the  pyramid  and  then  increasing  the  number  of  bins  used  by 
the  quantizer  for  each  level  of  the  pyramid  processed,  permitted  better  discrimination  and 
better  approximation  of  the  true  value.  Using  soft  thresholding  permits  the  storage  capacity 


of  the  token  to  be  maximized. 
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Use  of  the  binary  tree  storage  technique  from  our  model  provided  only  modest  color 
compression  on  the  order  of  15: 1  to  20: 1 .  More  importantly  though,  when  the  compression 
engine  was  replaced  with  a  model  that  used  a  more  sophisticated  two  dimensional  transform 
with  run  length  encoding,  it  was  found  that  the  color  separation  transforms  could  be  used  to 
achieve  average  compression  ratios  of  60: 1  with  only  minor  visual  defects.  It  was  possible 
to  reach  ratios  of  130: 1  with  only  modest  degradation  in  some  images  before  the  structural 
overhead  of  the  storage  method  limited  the  performance. 

A  suite  of  ten  varied  color  images,  all  512  x  512  and  24-bit  color,  were  used  for 
testing.  The  worst  compression  was  40: 1,  albiet  with  excellent  reconstruction.  The  image 
was  of  a  tiger  and  the  fur  of  the  tiger  acts  like  random  noise.  This  created  a  lot  of  detail 
coefficients  the  transformed  image  which  reduced  the  compressibility.  The  best  compression 
was  158:1  with  some  noticeable,  but  not  unreasonable  image  degradation.  This  image  was 
of  a  sunset  in  clouds.  Detail  in  the  clouds  was  lost  at  this  high  of  a  compression  ratio. 

We  also  tested  our  method  on  video  images.  Using  captured  video  frames, 
compression  ratios  of  70: 1  to  80: 1  were  easily  achieved.  The  video  images  suffer  from 
some  visual  defects  that  are  not  present  in  the  still  images.  These  defects  are  a  result  of 
motion  blur  and  of  the  poorer  resolution  of  video  signals  from  video  tape  and  tend  to  appear 
in  the  edge  detail  of  the  image.  Since  these  losses  complements  those  normally  generated 
by  the  wavelet  transform  we  were  able  to  get  excellent  results. 
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These  color  compression  ratios  easily  exceed  the  20: 1  compression  provided  by  the 
JPEG  standard.  With  further  research  on  the  optimal  color  and  wavelet  transforms,  and 
standardization  of  the  underlying  encoding  structure,  wavelet  transformation  provides  a 
viable  and  attractive  alternative  to  supercede  JPEG. 


Chapter  6 


Conclusion 

We  now  conclude  this  thesis  by  summarizing  the  major  points  and  providing  some 
directions  for  future  research. 


Summary 

In  the  thesis  we  have  presented  a  survey  of  the  major  forms  of  spatial  compression 
techniques  along  with  exploring  the  transform  domain  techniques.  We  have  seen  that  there 
is  a  variety  of  wavelet  transform  techniques  and  discussed  the  basic  operating  principles. 

We  presented  a  flexible  wavelet  transform  model  that  permits  compression  techniques 
to  be  explored  on  a  variety  of  image  sizes,  either  in  greyscale  or  color.  Although  the  binary 
tree  model  did  not  implement  the  most  effective  of  the  transform  methods  discussed,  it  does 
allow  the  compression  techniques  to  be  explored  on  a  range  of  images. 
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It  was  found  that  color  image  compression  can  be  effectively  implemented  by  using 
wavelet  transform  techniques  and  can  surpass  the  techniques  that  are  in  common  use  today 
in  the  JPEG  standard.  The  keys  to  successful  color  image  compression  are  transforming  the 
image  into  a  representation  that  minimizes  the  channel  bandwidth  required  to  reconstruct  the 
RGB  color  planes  of  the  image.  By  using  models  that  represent  the  luminance  and  color 
opposition  channels  effectively,  lossy  compression  can  be  tolerated  at  levels  beyond  lossy 
monochrome  compression  because  the  visual  degradation  is  masked  by  color  information. 

It  was  also  determined  that  another  key  to  successful  compression  is  the  selection  of 
effective  quantization  and  storage  techniques. 

Further  Research 

Although  the  binary  tree  storage  representation  implemented  was  not  as  effective  as 
is  expected,  it  is  known  that  this  is  in  part  due  to  the  inadequate  quantization  techniques.  If 
better  quantization  techniques  are  brought  to  bear  on  the  problem  and  heuristic  thresholding 
techniques  are  used,  it  is  believed  that  binary  tree  representation  of  the  transformed  image 
can  surpass  standard  run-length  encoding  of  the  image.  Improvements  in  the  quantization 
of  the  data  can  be  made  by  minimizing  the  mean  square  error  or  embedding  the  quantizer 
after  determining  the  significance  of  a  coefficient. 
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Another  viable  solution  that  was  not  explored  with  perhaps  better  discrimination  would 
be  to  first  soft  threshold  the  transformed  data,  build  the  binary  tree,  and  then  quantize  the 
remaining  data  based  on  the  only  the  retained  coeffients.  This  would  eliminate  from 
consideration  both  isolated  coefficients  and  coefficients  below  the  threshold  when 
constructing  the  quantization  bins. 

Human  visual  models  that  minimize  the  detectable  degradation  of  an  image  must  be 
further  explored.  Due  to  time  limitations,  the  logarithmic  transform  was  not  tested  against 
the  other  color  transfoms  to  determine  if  it  creates  a  more  visually  appealing  image  at  higher 
compression  ratios.  If  appropriate  transforms  are  developed  that  allow  quantization  errors 
to  be  greatest  within  ranges  that  are  less  detectable  by  the  human  eye  then  compression  ratios 
can  reach  even  greater  levels. 

Further  research  must  be  done  on  the  underlying  the  underlying  structure  of  the  storage 
method.  As  we  reach  greater  levels  of  compression,  the  actual  structure  overhead  then 
becomes  the  limiting  factor.  Some  of  the  structure  can  be  hard  coded  in  the  compression 
model  rather  than  being  passed  along  with  the  data.  Reducing  the  precision  of  the  control 
information  passed  along  in  the  model  from  32-bit  floating  point  values  to  16-bit  or  8-bit 
values  can  also  be  used  to  reduce  overhead. 

Lastly,  it  must  be  determined  how  much  loss  in  the  image  can  be  tolerated  with  real 
time  video.  The  movement  between  frames  and  the  30  frames/second  refreshing  of  a  video 
image  limits  the  eye’s  ability  to  detect  flaws,  it  may  be  possible  to  sustain  a  compression 
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ratio  of  over  100: 1  in  video  compression  without  significantly  affecting  the  visual  quality. 
Also,  methods  of  intraframe  compression  can  be  explored  by  only  saving  and  compressing 
the  differences  between  the  frames  instead  of  the  entire  next  frame. 

Implementation  of  the  wavelet  transform  in  hardware  should  perform  better  and  faster 
than  today’s  JPEG  because  of  the  reduced  complexity  of  the  wavelet  transform  over  the 
DCT.  By  using  a  pyramid  wavelet  transform  it  will  take  0(N)  versus  the  0(N  log  N)  of  the 


DCT. 


Appendix  I 


Computer  Program  Listings 


/* - 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  RGBKKK.C 

PURPOSE:  Splits  a  RGB  image  into  three  planes,  a  luminance 
plane  and  two  color  opposition  planes 

Input  file  is  a  RGB  (24-bit)  image  file 
3  output  files  are  created  each  is  a  float  file 


FEATURES:  can  handle  windows  BMP  format  files  also 


#include  <stdio.h> 
#include  <math.h> 
#include  <string.h> 
#include  <stdlib.h> 


^*******'fr**^*********|^^j^*****************^*******************^ 

int  main(int  argc,char  *argv[]) 

{ 

FILE  *in,  *R,  *G,  *B,  *Y,  *1,  *Q,  *YT,  *IT,  *QT: 
int  i,j,bmp’=0; 

char  filename[1 2],othername[1 2]={0,0,0,0,0,0,0,0,0,0,0,0}, 
name[8]={0,0,0,0,0,0,0,0}, 
ext[4]={0,0,0,0}, 
headerbuf[100]: 

unsigned  char  RGB[3],KKK2[3]: 
float  KKK[3]; 
int  width,height: 
unsigned  long  count=0: 
long  size; 
int  textout=0; 
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if  (argc<2)  printf("usage;  pgm  filename  [t]  {for  text  output}\n"); 
if  (argc>2&&argv[2][0]=='t')  textout=1 ; 
strcpy(name,argv[1]); 

if  ((in  =  fopen(argv[1],"rb"))  ==  NULL) 

{ 

fprintf(stderr, "Input  %s  not  foundAn", filename): 
exit(1): 

} 

for  (i=0;i<=14;i++)  othername[i]=0; 

strcpy(othername,argv[1]); 

strncat(othername,".res",4): 

if  (strcmp(ext,'’bmp")==0)  /*  assume  BMP  file  */ 

{ 

fread(headerbuf,1 ,57, in); 

if  (headerbuf[0]!='B'll  headerbuf[1]!='M') 

{ 

printf("File  not  windows  bitmap\n"); 
exit(O); 

} 

if  (headerbuf[28]!=24) 

{ 

printf("lmage  is  not  24  bit  color  windows  bitmap\n"); 
exit(O); 

} 


bmp  =1; 

width  =  (((int)  headerbuf[19])«8)  +  (int)  headerbuf[1 8]; 
height  =  (((int)  headerbuf[23])«8)  +  (int)  headerbuf[22]; 

} 

else  /*  assume  RGB  file  */ 

{ 

size=filesize(in): 

height=wldth=(int)sqrt((double)size/3); 

if((long)height*width*3-size!=0) 

{ 

printf("filesize  error\n“); 
exit(1); 

} 

} 

printf("height=%i,  width=%i,  size=%lu\n",height,width,(unsigned  long)height*width); 

for  (i=0;i<=14:i++)  othername[i]=0; 

strcpy(othername,name): 

strncat(othername,".k1f“,4); 

if  ((Y  =  fopen(othername,“wb"))  ==  NULL) 

{ 

fprintf(stderr, "Output  file  %s  not  createdAn",othername); 
exit(1); 
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for  (i=0;i<=14:i++)  othername[i]=0: 

strcpy(othername,name); 

strncat(othername,''.k2f",4): 

if  ((I  =  fopen(othername,''wb“))  ==  NULL) 

{ 

fprintf(stderr, "Output  file  %s  not  created.\n'',othername): 
exit(1): 

} 

for  (i=0;i<=14:i++)  othername[i]=0; 

strcpy(othername,name); 

strncat(othername,".k3f'',4); 

if  ({Q  =  fopen(othername,''wb"))  ==  NULL) 

{ 

fprintf(stderr, "Output  file  %s  not  created.\n",othernanne); 
exit(1); 

} 

for(i=0;i<height*width:i++) 

{ 

fread(RGB,1,3,in): 

if  (bmp) 

{ 

KKK[0]=(0.575*(float)RGB[2]+0.615*(float)RGB[1]+0.540*(float)RGB[0]): 

KKK[1]=(0.608‘(float)RGB[2]+0.120*(float)RGB[1]-0.785*{float)RGB[0]); 

KKK[2]=(0.548*(float)RGB[2]-0.779‘(float)RGB[1]+0.305*(float)RGB[0]); 

} 

else 

{ 

KKK[0]=(0.575‘(float)RGB[0]+0.615*(float)RGB[1]+0.540*(float)RGB[2]): 

KKK[1]=(0.608*(float)RGB[0]+0.120*(float)RGB[1]-0.785*(float)RGB[2]); 

KKK[2]=(0.548*(float)RGB[0]-0.779*(float)RGB[1]+0.305*(float)RGB[2]): 

} 

/*  this  is  for  char  based  files,  scaling  and  tranlation  functions  to  fit  back  into  char  files 
KKK2[0]={unsigned  char)(KKK[0]/1.73): 

KKK2[1]=(unsigned  char)(KKK[1]/1.513+132.3): 

KKK2[2]=(unsigned  char)(KKK[2]/1 .632+1 21.7); 
fwrite(&KKK2I0],1,1,Y); 
fwrite(&KKK2[1], 1,1,1); 
fwrite(&KKK2[2],1,1,Q); 

7 

fwrite{&KKK[0]  ,sizeof  (float) ,  1 ,  Y) ; 
fwrite(&KKK[1],sizeof(float),1,l); 
fwrite(&KKK[2] ,  sizeof  (float)  ,1,0); 

if  (textout) 

{ 

fprintf(YT,"%6.2f  ",KKK[0]); 


fprintf(IT,"%6.2f  ",KKK[1]): 
fprintf(QT,''%6.2f  '',KKK[2]): 

} 

if  (count++%width==0) 

{ 

printf(’'%lu\r", count): 

} 

} 

printf("\n’'); 

fclose(in); 

fclose(Y): 

fclose(l): 

fclose{Q): 

return  0; 

} 


/* . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  KKKRGB.C 

PURPOSE:  Reconstructs  an  RGB  file  from  3  float  files,  recovers 
RGBKKK’s  transform 


#include  <stdio.h> 
#include  <math.h> 
#include  <string.h> 
#include  <stdlib.h> 


^********************|^^j^***********^************,^********^^,j^* ! 

int  main(int  argc,char  *argv[]) 

{ 

FILE  *out,  *K1,  *K2,  *K3: 
int  i,j; 

char  filename[20],othername[12]={0,0,0,0,0,0,0,0,0,0,0,0}, 
name[8]={0,0,0,0,0,0,0,0}, 
ext[4]={0,0,0,0}, 
headerbuf[100]; 
unsigned  char  KKK2[3]; 
float  RGB[3],KKKt3]; 
int  width,height: 
unsigned  long  count=0; 


if  ((out  =  fopen(argv[1],"wb"))  ==  NULL) 

{ 

fprintf(stderr, "output  %s  not  created.\n'',filename); 
exit(1); 

} 

for  (i=0;i<=14;i++)  othername[i]=0; 
strcpy(othername,argv[1  ]); 
strncat(othername,".k1f'',4): 

if  ((K1  =  fopen(othername,''rb”))  ==  NULL) 

{ 

fprintf(stderr, "input  file  %s  not  found.\n“,othername); 
exit(1); 

} 

for  (i=0;i<=14;i++)  othername[i]=0; 
strcpy(othername,argv[1  ]); 
strncat(othername,“.k2f",4); 

if  ((K2  =  fopen(othername,''rb“))  ==  NULL) 

{ 
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fprintf(stderr, "input  file  %s  not  found.\n",othername); 
exit(1): 

} 

for  (i=0:i<=14:i++)  othername[i]=0; 

strcpy(othername,argv[1]): 

strncat(othername,".k3f",4): 

if  ((K3  =  fopen(othername,''rb“))  ==  NULL) 

{ 

fprintf(stderr, "input  file  %s  not  found.\n'',othername); 
exit(1): 

} 

while  (!feof(K1)) 

{ 

fread(&RGB[0],sizeof(float),1  ,K1 ); 
fread(&RGB[1],sizeof(float),1,K2): 
fread(&RGB[2],sizeof(float),1,K3); 
if (feof(K1 ))  break; 


KKK[0]=0.575*RGB[0]+0.608*RGB[1]+0.547*RGB[2]; 
KKK[1  ]=0.61 5*RGB[0]+0. 1 20*RGB[1  ]-0.779‘RGB[2]; 
KKK[2]=0.539*RGB[0]-0.785*RGBt1]+0.305*RGB[2]; 

KKK2[0]=(unsigned  char)(KKK[0]+.5); 

KKK2[1  ]=(unsigned  char)(KKK[1  ]+.5); 
KKK2[2]=(unsigned  char)(KKK[2]+.5); 

fwrite(KKK2,1,3,out): 

if  (count++%1000==0) 

{ 

printf("%lu\r", count); 

} 

} 

printf("\n"); 
fclose(out); 
fclose(K1); 
fclose(K2); 
fclose(K3); 
return  0; 

} 


/* . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  ZT2.C 

PURPOSE:  Builds/recovers  a  zero  subtree  from  a  wavelet 
transform  file 

Input  file  is  a  binary  float  file 
output  file  is  a  char  file 

Expects  a  data  header  of  1  +  2*levels  processed  (floats) 

Containing  levels  processed  and  mean  and  standard  deviation  of  each 
level  processed. 

FEATURES:  intelligent  selection  of  data  nodes  to  keep  uses 
a  2-bit  binary  prefix  to  encode  directional  decision  data 

CONSTRAINTS:  no  file  size  detection  routines 


#include  <stdio.h> 
#include  <stdlib.h> 
#include  <math.h> 
#include  <string.h> 

#define  DEBUG  0 


GLOBALS 

★***W**'**************w***iir****ifc*^******w*******************^ 

FILE  *infile,  *outfile; 
int  outcount  =  0; 
int  count!  26  =  0; 
int  count!  27  =  0; 
int  countm!27=  0; 
int  countspec=  0; 
int  countlong=  0; 
int  threshold!  =0; 
float  threshold2  =  0; 
int  lookahead  =  2; 
int  hsize  =  -!; 
int  level  =  0; 
int  direction  =  ! ; 
int  info  =  ! ; 
int  dbgcount=!00: 
int  dbg=0: 
int  DEBUG2  =  0; 


void  setoptions(int  argc.char  *argv[]) 

{ 


if(argc  <  4) 

{ 

printf(”Zero  T ree  -  zero  tree  for  a  wavelet  transform  file  \n"); 
printfC'This  program  is  for  float  format  files  \n"); 
printfC'Usage;  pgm  filename  LEVELS  DIR  [options]\n''); 
printf(“DIR  1  forward  -1  backward  \n"); 
printf(“LEVELS  levels  to  process  \n"); 
printf("options:  \n"); 

printf(''-T#  primary  threshold  {  0  default }  \n"); 
printf("-t#  secondary  threshold  { 0  default }  \n"); 
printf("-a#  lookahead  levels  {  2  default  }\n''); 
printf("-h#  size  of  header  { checks  file  by  default }  \n"): 
exit(O): 

} 

for  (i=4;i<argc;i++) 

{ 

switch  (argv[i][1]) 

{ 

case  'a'rlookahead  =  atoi(&argv[i][2]): 
break; 

case  T':threshold1  =  atoi(&argv[i][2]); 
break; 

case  't':threshold2  =  atof(&argv[i][2]); 
break; 

case  'h'rhsize  =  atoi(&argv[i][2]); 
info  =  0; 
break; 

case  'd';DEBUG2  =  atoi(&argv[i][2]); 
break; 

} 

} 

level  =  atoi(argv[2]); 

direction  =  atoi(argv[3]); 

if  (!abs(direction)==1)  direction=  1; 

if  (direction!=1)  printf("levels  =  %i,  direction  %i\n'', level, direction); 
else  printf("levels  =  %i,  direction  %i,  lookahead  =  %i,  threshi  =  %\," 
“  thresh2  =  %i\n“, level, direction, lookahead, thresholdl, 
threshold2); 


void  *  getvect(int  n,int  size) 

{ 


int  i; 

void  *vectptr; 


vectptr  =  (void  *)  calloc(n,size); 
if(vectptr  ==  NULL)  { 
printfC'Error:  dynamic  memory  fail\n"); 
exit(O); 


} 

return(vectptr); 

} 


^'Ar********************:^***************************<r*5lf***** 

r  file  size  determinator  7 
int  filesize(FILE  ‘stream) 

{ 

int  curpos,  length; 
curpos  =  ftell(stream); 
fseek(stream,  OL,  SEEK.END); 
length  =  ftell(stream); 
fseek(stream,  curpos,  SEEK_SET): 
return  length; 

} 


signed  char  *encode(signed  char  ‘treeptr, float  ‘input) 

{ 


int  lowbyte=Oxff; 


if{abs({int)‘input)<128)  /‘  -128  is  reserved  code‘/ 

{ 

‘treeptr+H-  =  (signed  char)‘input; 
outcount++; 

} 

else  if  (abs((int)‘input)<=32895) 

{ 

‘treeptr++=-128; 

‘treeptr++=  (signed  char)((int)‘input»8); 

‘treeptr++=  (unsigned  char)((int)‘input&lowbyte); 

outcount+=3; 

countlong++; 

if(DEBUG)  printf("decomp  %i  into  %i  +  %i\n",(int)‘input, 
(int)‘(treeptr-2),(int)‘(treeptr-1)); 

} 

else 

{ 

printfC'ABORTED:  range  error  exceeded  on  abs  calc  =  %i\n",(int)‘input); 
fwrite(&treeptr[-outcount],sizeof(signed  char),outcount,outfile); 
printf("%i  token  saved  in  outfile\n'',outcount); 
exit(O); 

} 

return  treeptr; 

} 


int  threshcheck(float  ‘loc[],int  level, int  depth,int  retval) 

{ 

int  i,  chkval, breadth; 
breadth  =  1«(depth-1); 


chkval  =  (threshold1-(int)(threshold2*depth>0))  ? 
threshold1+(int)(threshold2*depth):0; 

if  (!(retval&0x02))  /*  don't  repeat  check  if  already  found  one  */ 

{ 

for  (i=0;i<breadth;i++)  /‘left  side*/ 

{ 

if{abs(loc[level][i])>=chkval) 

{ 

retvall=0x02: 

if  (DEBUG  )  printfC'left  trigger  at  lev  %i,  off  %i,  ckval  %l," 

"  val  %i\n", level,!, chkval, (int)loc[level][i]); 

} 

} 

} 

if  (!(retval&0x01))  /*  don't  repeat  check  if  already  found  one  */ 

{ 

for  (i=breadth;i<breadth*2;i++)  /*  right  side*/ 

{ 

if(abs(loc[level][i])>=chkval) 

{ 

retval  1=0x01 ; 

if  (DEBUG  )  printfC'right  trigger  at  lev  %i,  off  %i,  ckval  %i," 
"  val  %i\n", level, l,chkval,(int)loc[level][i]); 

} 

} 

} 

if  (retval==3)  return  retval;  /*  both  sides  signif,  stop  */ 
else  if  (depth<lookahead&&level>1 ) 

return  threshcheck(loc,level-1  ,depth+1  .retval); 
else  return  retval; 

} 


signed  char  *treetraverse(float  *vector,signed  char  ‘treeptr,  float  ‘locQ, 
float  *tree[],int  level) 

{ 

int  i; 

int  nextievel; 

Int  initval  =0; 
int  initdepth  =1 ; 
signed  char  temp; 


/‘encode  left  node  and  subtree  if  not  zero*/ 


if  (/*abs(*loc[ievel])<=threshold1  &&*/level>1 ) 

{  /*  need  to  see  if  anything  sig  is  below*/ 

switch(threshcheck(loc, level-1  .initdepth, initval)) 

{ 

case  0:  /*  both  subtrees  insignificant  */ 
temp  =  (signed  char)  *loc[level];  /*  get  input  */ 


if  (abs(*loc[level])>31 )  /*  if  beyond  range,  special  code  */ 


{ 


*treeptr++  =  OxeO;  /*  dir  bits  and  output  range  exception  */ 

temp  =  treeptr[-1  ];  /*  looks  like  1 1 1 00000  V 

*treeptr++  =  (signed  char)  *loc[level]; 

countspec++: 

outcount+=2: 

} 

else  /*  fits  in  range  */ 

{  /*  1 00000000  ->  001 00000  */ 
if  (temp&0x80)  temp  l=  0x20;  /*  move  the  sign  bit  */ 

temp  l=  OxcO;  /*  mask  directional  bits  11000000  */ 
*treeptr++  =  temp;  /*  save  encoded  char  */ 

outcount++; 

} 

count126++: 

if  (DEBUG2->0) 

{ 

printfC'level  %i,  stop  at  %i,  i=%i,  o=%i  -  both  pruned, " 
"going  back\n“,  level, (int)(loc[level]-tree[level]), 
(int)(signed  char)*loc[level],(int)temp); 


nextievel  =  level;  /‘update  subtree  pointers  2  ea  level*/ 
for{i=nextlevel;i>0:i--) 

{ 

loc[i]  +=  1«(nextlevel-i); 

} 

return  treeptr; 

case  1 :  /*  left  subtree  insignificant,  process  right*/ 
temp  =  (signed  char)  *loc[level];  /*  get  input  */ 

if  (abs(*loc[level])>31)  /*  if  beyond  range,  special  code  */ 

{ 

*treeptr++  =  OxaO;  /*  dir  bits  and  output  range  exception  */ 

temp  =  treeptr[-1]:  /*  looks  like  10100000  */ 

*treeptr++  =  (signed  char)  *loc[level]; 

countspec++; 

outcount+=2; 

}  . 

else  /*  fits  in  range  */ 

{  /*  1 00000000  ->001 00000  */ 

if  (temp&0x80)  temp  1=  0x20;  /*  move  the  sign  bit  */ 

temp  &=  0x3f;  /*  clear  directional  bits  001 11111  */ 

temp  1=  0x80;  /*  mask  directional  bits  10000000  */ 

*treeptr+-i-  =  temp;  /*  save  encoded  char  */ 

outcount-n-; 

} 

count127-(-i-; 

if  (DEBUG2->0) 

{ 

printfC'level  %i,  right  at  %i,  i=%i,  o=%i  -  pruned  left, " 


} 


'going  right  \n",  level, (int)(loc[level]-tree[level]), 
(int)(signed  char)*loc[level],(int)temp): 


/‘update  subtree  pointers,  before  doing  right  7 
loc[level]++;  /*  consume  input  7 

nextievel  =  ievel-1 ; 
for(i=nextlevel;i>0;i--) 

{ 

loc[i]  +=  1«(nextlevel-i); 

} 

treeptr  =  treetraverse(vector,treeptr,loc,tree,level-1); 
return  treeptr; 

case  2;  /*  right  subtree  insignificant,  process  Ieft7 
temp  =  (signed  char)  *loc[level];  /*  get  input  7 

if  (abs(*loc[level])>31)  /*  if  beyond  range,  special  code  7 

{ 

*treeptr++  =  0x60;  /*  dir  bits  and  output  range  exception  7 

temp  =  treeptr[-1  ];  /*  looks  like  01 1 00000  7 

*treeptr++  =  (signed  char)  *loc[level]; 

countspec++; 

outcount+=2; 

} 

else  r  fits  in  range  7 

{  /*  1 00000000  ->  001 00000  7 

if  (temp&0x80)  temp  1=  0x20;  /*  move  the  sign  bit  7 

temp  &=  0x3f;  /*  clear  directional  bits  001 11111  7 

temp  1=  0x40;  /*  mask  directional  bits  01000000  7 

*treeptr++  =  temp;  /*  save  encoded  char  7 

outcount++; 

} 

if  (DEBUG2->0) 

{ 

printf("level  %i,  left  at  %i,  i=%i,  o=%i  -  pruned  right, " 

"going  left  \n",  level, (int)(loc[level]-tree[level]), 
(int)(signed  char)*loc[level],(int)temp); 

} 

countm127++; 

treeptr  =  treetraverse(vector,treeptr,loc,tree,level-1); 

if  (DEBUG2-->0)  printf("level  %i,  right  at  %i  skipped\n'',level, 
(int)(loc[level]-tree[level])); 

/‘update  subtree  pointers,  after  left  done  ‘/ 
nextievel  =  level- 1 ; 
f  or(i=nextlevel;  i>0;  i-) 

{ 

loc[i]  +=  1«(nextlevel-i); 

} 

loc[level]++; 
return  treeptr; 


/‘  case  3;  both  sides  are  significant  fall  through  and  process  ‘/ 
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} 


} 


r  or  else  process  normally  V 

temp  =  (signed  char)  *loc[level]; 

if  (DEBUG>0)  printfC'in  temp  =  %i,",(int)temp); 

if  (abs(*loc[level])>31)  /*  if  beyond  range,  speciai  code  */ 

{ 

*treeptr++  =  0x20;  /*  dir  bits  and  output  range  exception  7 

temp  =  treeptr[-1  ];  /*  looks  like  001 00000  7 

if  (DEBUG>0)  printfC  out  code  =  %i,  '',(int)temp); 

*treeptr++  =  (signed  char)  *loc[level]: 

countspec++; 

outcount+=2; 

if  (DEBUG>0)  printf("out  temp  =  %i\n",(int)temp); 

} 

else  r  fits  in  range  7 

{  r  1 00000000  ->  001 OOOOO  7 

if  (temp&0x80)  temp  1=  0x20;  /*  move  the  sign  bit  7 

temp  &=  0x3f;  /*  clear  directional  bits  001 11111  7 

*treeptr++  =  temp;  /*  save  encoded  char  7 

outcount++; 

if  (DEBUG>0)  printf(“  out  temp  =  %i  FITS\n",(int)temp): 

} 

if(level>1 ) 

{ 


if  (DEBUG2-->0)  printf(''level  %i,  left  at  %i,  i=%i,  o=%i\n", 
level, loc[level]-tree[level],(int)((signed  char)*loc[level]), 
(int)((signed  char)*(treeptr-1 ))); 

treeptr  =  treetraverse(vector,treeptr,loc,tree,level-1); 

/‘encode  right  subtree  if  level>1 7 

if  (DEBUG2-->0)  printf("level  %i,  right  at  %i\n“,level, 
loc[level]-tree[level]); 

treeptr  =  treetraverse(vector,treeptr,loc,tree,level-1); 

} 

else 

{. 

if  (DEBUG2->0)  printf("level  %i,  done  at  %i,  i=%i,  o=%i\n", 
level, loc[level]-tree[level],(int)((signed  char)*loc[level]), 
(int)((signed  char)*(treeptr-1 ))); 

} 

loc[level]++;  /*  update  location  7 
return  treeptr; 

} 


^★***************:*;***!fc*****************^*****^f*Hr******Ht<r****************^ 


void  zerotree(float  ‘vector, signed  char  *treevector,int  size.int  levels) 

{ 

int  i,  maxlevel; 
float  **loc; 
float  “tree; 
signed  char  ‘treeptr; 

printf("setting  up  tree  pointersNn"); 

i=0; 

while(1«i<size)  maxlevel=++i; 
loc  =  getvect(maxlevel+1,sizeof(int)); 
tree  =  getvect(maxlevel+1  ,sizeof(int)); 

printf ("maxlevel  is  %i,  levels  requested  are  %i\n'', maxlevel, levels); 

if(levels>=maxlevel) 

{ 

printf  ("too  many  levels  requestedNn"); 
exit(O); 

} 

tree[maxlevel]  =  vector;  /*  should  add  test  for  1>  2**levels  7 
loc[maxlevel]  =  vector; 

for  (i=0;i<maxlevel;i++)  /*  set  up  pointers  0  =  end,  maxlevel  =  startV 

{ 

tree[i]  =  &vector[size»i] ; 
loc[i]  =  tree[i]; 

} 

if(DEBUG)printf("vector  address  =  %i\n",(int)vector); 

for  (i=maxlevel;i>0;i--) 

{ 

printf("  loc  %i  offset  %i  contains  %i  floats", 
i,loc[i]-loc[maxlevel],loc[i-1]-loc[i]); 
if(DEBUG)  printfC,  starts  at  addr%i",(int)loc[i]); 
printf("\n"); 

} 

treeptr  =  treevector; 
i=0; 

while  (&vector[i]<  tree[levels]) 

{ 

treeptr  =  encode(treeptr,&vector[i]); 
if(DEBUG)  printf("vect  =  %6.2f  gives  c  =%i  \n", 
vector[i],(int)treeptr{-1  ]); 
i++: 

} 

printf("%i  tokens  copied\n",i); 

printf ("%i  tokens  to  process\n",tree(levels-1]-loc[levels]); 


while  (loc[leveis]<tree[levels-1])  /*  begin  recursive  tree  search  call  */ 

{ 

treeptr  =  treetraverse(vector,treeptr,loc, tree, levels); 

} 


signed  char  *decode(signed  char  ‘treeptr, float  ‘outputloc) 

{ 

int  temp; 

if  (*treeptr==-128)  /*  rebuild  integer  value  */ 

{ 

treeptr++; 

temp  =  (int)((signed  char)*treeptr++); 
temp  «=8; 

temp  +=  (int)((unsigned  char)*treeptr++); 

‘outputloc  =  (float)temp; 
outcount++; 

if  (DEBUG)  printfC'rebuild  %i  from  %i  +  %i\n", 

(int)*outputloc,(int)*(treeptr-2),(int)*(treeptr-1)); 

} 

else 

{ 

*outputloc=(float)*treeptr++; 

outcount++; 

} 

return  treeptr; 

} 


y***************************************y|.********1lr******************'****  j 

signed  char  *buildtree(float  *tree[],  signed  char  *treevector,int  size, 
signed  char  ‘treeptr,  int  level, int  offset) 

{ 

signed  char  ‘treelast,  ctemp,  direction; 
int  temp; 

if  (offset>=tree£ievel-1]-tree[level])  /‘  if  input  exhasted  go  back  ‘/ 

{ 

printfC'level  %i  ,si2e  %i,  index  %i  completed\n'', 
level,tree[level-1]-tree[level],offset); 
return  treeptr; 

} 

r  get  direction  and  then  decode  char  */ 
direction  =  (‘treeptr»6)&0x03;  /‘  shift  6,  mask  with  0000001 1  ‘/ 

ctemp  =  ‘(treeptr++)&0x3f;  /‘  mask  out  direction  001 11111  ‘/ 
outcount++;  /*  consume  token  ‘/ 

if(DEBUG>0)  printfC'ctempin  =  %i,  ",(int)ctemp); 
if  (ctemp==0x20)  r  then  special  code  exception,  rebuild  number  ‘/ 

{ 

tree[levei][offset]  =  ‘treeptr++;  /‘  output  whole  next  char  ‘/ 

/‘  consume  token  ‘/ 


} 


else  r  convert  this  code  back  */ 

{ 

if(ctemp&0x20)  ctemp  =  ctemplOxeO;  /*  if  expand  sign  bit  */ 
tree[level][offset]  =  (float)  ctemp;  /*  output  value  of  token  7 

} 

if  (DEBUG>0)  printf{"*treeptr=  %i,  direction  =  %i,  ctemp  =  %i\n", 
{int)*(treeptr-1 ),  (int)  direction,  (int)ctemp); 

switch  (direction)  /*  now  decide  which  direction  to  go  7 

{ 

case  3:  /*  1 1 XXXXXX  ->  0000001 1  stop  code  7 

if  (DEBUG2->0) 

{ 

printf("level  %i,  stop  at  %i,  i=%i,  o=%i  -  both  pruned," 

"  going  back  \n“, 

level, (int)offset,(int)(signed  char)'*(treeptr-1 ), 
(int)tree[level][offset]); 

} 

return  treeptr; 


case  2:  /*  1 0XXXXXX  ->  0000001 0  go  right  code,  skip  left  7 

if  (DEBUG2->0) 

{ 

printf("level  %i,  right  at  %i,  i=%i,  o=%i  -  pruned  left," 

"  going  right  \n", 

ievel,(int)offset,(int)(signed  char)*(treeptr-1 ), 
(int)tree[level][offset]); 

} 

if(level>1)  /‘continue  7 

{ 

treeptr  =  buildtree(tree,treevector,size,treeptr, 
level-1, 2*(offset)+1); 

} 

return  treeptr; 


case  1 :  .  /*  01  XXXXXX  ->  00000001  go  left  code,  skip  right  7 

if  (DEBUG2->0) 

{ 

printf("level  %i,  left  at  %i,  i=%i,  o=%i  -  prune  right," 

"  going  left  \n", 

level,(int)offset,(int)(signed  char)*(treeptr-1 ), 
(int)tree[level][offset]); 

} 

/‘go  left  if  level>1  ‘/ 
if(level>1) 

{ 

treeptr  =  buildtree(tree,treevector,size,treeptr,level-1, 
2*(offset)); 
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} 

if  (DEBUG2-->0)  printf("level  %i,  right  at  %i  skipped\n'', level, 
offset): 

return  treeptr; 


default:  /*  OOXXXXXX  ->  00000000  no  code,  process  normally  */ 

/*go  left  if  level>1 )  */ 

if(level>1) 

{ 


if  (DEBUG2-->0)  printf("level  %i,  left  at  %i,  i=%i,  o=%i\n", 
level,offset,(int)((signed  cha r) ‘(treeptr- 1)), 
(int)((float)tree[level][offset])): 

treeptr  =  buildtree(tree,treevector,size,treeptr, 
level-1 ,2*(offset)): 

/*  and  then  do  the  right  side  */ 


if  (DEBUG2->0)  printf("level  %i,  right  at  %i\n", 
level,offset); 

treeptr  =  buildtree(tree,treevector,size,treeptr, 
level-1 ,2*(offset)-h1 ); 

} 

else 

{ 

if  (DEBUG2->0)  printf("level  %i,  done  at  %i,  i=%i,  o=%i\n", 
level,offset,(int)((signed  char)*(treeptr-1)), 
(int)((float)tree[level][offset])): 

} 

} 

return  treeptr; 

} 


void  unzerotree(float  *vector,signed  char  *treevector,int  size,int  levels) 

{ 


int  i,  maxlevel; 
float  “loc; 
float  “tree; 
signed  char  ‘treeptr; 
int  offset  =  0; 


printfC'rebuilding  tree  pointers\n''); 
i=0: 

while(1«i<size)  maxlevel=-i-+i; 
loc  =  getvect(maxlevel-(-1  ,sizeof(int)); 


tree  =  getvect(maxlevel+1,sizeof(int)): 

printf("maxlevel  is  %i,  levels  to  restore  are  %i\n'',maxlevel, levels); 


if(levels>=maxlevel) 

{ 

printfC'too  many  levels  requested\n"): 
exit(O); 

} 

tree[maxlevel]  =  vector; 
loc[maxlevel]  =  vector; 

for  (i=0;i<maxlevel;i++)  /‘set  up  pointers  0  =  end,  level+1  =  begining  */  { 
tree[i]  =  &vector[size»i] ; 
loc[i]  =  tree[i]; 

} 

if  (DEBUG)  printf("vector  address  =  %i\n",(int)vector): 
for  (i=maxlevel;i>0;i--) 

{ 

printf("  loc  %i  offset  %i  contains  %i  floats", 
i,loc[i]-loc[maxlevel],loc[i-1]-loc[i]); 
if  (DEBUG)  printf(",  starts  at  addr  %i",(int)  loc[i]); 
printf(''\n"); 

} 

treeptr  =  treevector; 
i=0; 

while  (&vector[i]<  tree[levels]) 

{ 

treeptr  =  decode(treeptr,&vector[i]); 
if(DEBUG)  printf("vect  =  %6.2f  \n",vector[i]); 
i++; 

} 

printf("%i  tokens  copied\n",i); 
offset=0; 

/*  while  there  is  still  input*/ 
while(offset<=tree[levels]-tree[maxlevel]) 

{ 

treeptr  =  buildtree(tree,treevector,size,treeptr,levels,offset++); 

} 

} 


^★★★**************************^*i(f***************ikr***************^Hk******^ 

void  main(int  argc,  char  *argv[]) 

{ 

int  fsize; 

int  xdim=0,  ydim=0,  size=0; 
float  ‘vector; 


float  *headvector; 
signed  char  ‘treevector; 
float  flevelq,  flevelzt; 

char  string[25]={0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0}; 

setoptions(argc,argv); 
infile  =  fopen(argv[1],"rb"); 
if(infile  ==  NULL) 

{ 

printfC'File  open  failure\n"); 
exit(O); 

} 

strcpy(string,argv[1]): 
if  (direction>0)  strcat(string,".zt2"); 
else  strcat(string,".tz2"); 
outfile  =  fopen(string,''wb''); 
if(outfile==  NULL) 

{ 

printfC'File  open  failure\n"): 
exit(O); 

} 

size  =  filesize(infile); 
if  (hsize==-1)  /*  no  header  specified  V 
{ 

fread(&flevelq,sizeof(float),1, infile):  /‘levels  processed  by  quant*/ 
hsize  =  (flevelq+info)*2; 

} 

if(direction>0)  size  =  ((size-1  )/4)-hsize; 
else  size  -=  (hsize+info); 

headvector=getvect(hsize,sizeof(float)); 
treevector  =getvect(size*2,sizeof(signed  char)); 

if  (direction==1) 

{ 

vector  =  getvect(size,sizeof(float)); 

fwrite(&flevelq,sizeof(float),info,outfile);/*levels  processed  by  quant*/ 
flevelzt  =  (float)level; 

fwrite(&flevelzt,sizeof(float),1,outfile);/*levels  process  by  zt  */ 
fread(headvector,sizeof(float),hsize, infile);  /*get  data  vector*/ 
fwrite(&size,sizeof(float),1, outfile): 
printfC'saving  header  of  size  %i\n",hsize); 
fwrite(headvector,sizeof(float),hsize,outfile);  /‘save  data  vector*/ 
fread(vector,sizeof(float),size,infile):  /‘get  image  vector  */ 

zerotree(vector,treevector,size,level); 
printf("output  file  write  of  %i\n", 
fwrite(treevector,sizeof(signed  char),outcount,outfile)); 
printf("%s  =  %i  tokens,  converted  to  %i  char\n",argv[1],size,outcount); 
printf(‘'%i  long  codes  output\n“,countlong); 

printfC'output  stop  =  %i,  left  =  %i,  right  =  %i,  spec  =  %i,  total  code  tokens  %i\n", 

counti  26,count1 27,countm1 27,countspec,count1 26-HCOuntl  27-i-countml  27-i-countspec); 
printf("token  overhead  %3.2f%7o\n". 
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1 00.0*{(f  loat)(count1 26+count1 27+countm  1 27+countspec)/(float)outcount)); 

else 

{ 

fread(&flevelzt,sizeof(float),1,infile);/*level  processed  by  zt  7 
level  =  (int)flevelzt: 
fread(&fsize,sizeof(float),1, infile): 
vector  =  getvect(fsize,sizeof(float)); 
fwrite(&flevelq,sizeof(float),info,outfile); 
fread(headvector,sizeof(fioat),hsize, infile);  /‘get  data  vectorV 
printfC'saving  header  of  size  %i\n“,hsize): 
fwrite(headvector,sizeof(float),hsize,outfile);  /‘save  data  vectorV 
fread(treevector,sizeof(signed  char), size, infile); 
printfC'reading  in  %i  data  tokens  \n“,size); 
unzerotree(vector,treevector,fsize, level); 

printf(''%s  =  %i  tokens,  converted  to  %i  floats\n",argv[1],size,outcount); 
printfC'output  file  write  of  %i\n", 
fwrite(vector,sizeof(float),fsize,outfile)): 

} 

free(vector); 

free(treevector); 

fclose(infile); 

fclose(outfile): 

} 
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/* . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  ZT.C 

PURPOSE:  Builds/recovers  a  zero  subtree  from  a  wavelet 
transform  file 

Input  file  is  a  binary  float  file 
output  file  is  a  char  file 

Expects  a  data  header  of  1  +  2*levels  processed  (floats) 

FEATURES:  intelligent  selection  of  data  nodes  to  keep,  uses  8  bit  stop  tolens 
CONSTRAINTS:  no  file  size  detection  routines 


#include  <stdio.h> 
#include  <stdlib.h> 
#include  <math.h> 
#include  <string.h> 

#define  DEBUG  0 
#define  DEBUG2  0 


^******************«»**********'*****************************Hf 

GLOBALS 

*******************^t***************************************^ 

FILE  *infile,  ‘outfile; 
int  outcount  =  0; 
intcount126  =  0: 
int  count127  =  0; 
int  countm127=  0; 
intthresholdl  =  0; 
int  threshold2  =  0; 
int  lookahead  =  2; 
int  hsize  =  -1 ; 
int  level  =  0; 
int  direction  =  1 ; 
int  info  =  1; 
int  dbgcount=100: 
int  dbg=0; 


void  setoptions(int  argc,char  *argv[]) 


{ 


int  i; 


if(argc  <  4) 

{ 

printfC'Zero  Tree  -  zero  tree  for  a  wavelet  transform  file  \n"); 


printfC'This  program  is  for  float  format  files  \n"); 
printf("Usage:  pgm  filename  LEVELS  DIR  [options]\n"); 
printf("DIR  1  forward  -1  backward  \n"); 
printf("LEVELS  levels  to  process  \n"); 
printf(“options:  \n''); 

printf("-T#  primary  threshold  {  0  default }  \n''); 
printf("-t#  secondary  threshold  {  0  default }  \n"); 
printf(''-a#  lookahead  levels  {  2  default  }\n''); 
printf(“-h#  size  of  header  {  checks  file  by  default }  \n''); 
exit(O); 

} 

for  (i=4;i<argc;i++) 

{ 

switch  (argv[i][1]) 

{ 

case  'a'rlookahead  =  atoi(&argv[i][2]): 
break; 

case  T':threshold1  =  atoi(&argv[i][2]); 
break; 

case  't':threshold2  =  atoi(&argv[i][2]); 
break; 

case  'h';hsize  =  atoi(&argv[i][2]); 
info  =  0; 
break; 

} 

} 

level  =  atoi(argv[2]); 

direction  =  atoi(argv[3]); 

if  (!abs(direction)==1)  direction=  1; 

if  (direction!=1)  printf(''levels  =  %i,  direction  %i\n'', level, direction); 
else  printfC'levels  =  %i,  direction  %i,  lookahead  =  %i,  threshf  =  %i," 
"  thresh2  =  %i\n'', level, direction, lookahead, thresholdl, 
threshold2); 


void  *  getvect(int  n,int  size) 

{ 

int  i; 

void  *vectptr; 


vectptr  =  (void  *)  calloc(n,size); 
if(vectptr  ==  NULL)  { 
printf("Error:  dynamic  memory  fail\n"); 
exit(O); 

} 

return  (vectptr); 

} 


/*  file  size  determinator  7 
int  filesize(FILE  ‘stream) 


{ 

int  curpos,  length; 
curpos  =  ftell(stream); 
fseek(stream,  OL,  SEEK_END); 
length  =  ftell(stream): 
fseek(stream,  curpos,  SEEK_SET); 
return  length; 

} 


signed  char  *encode(signed  char  ‘treeptr, float  ‘input) 

{ 

int  lowbyte=Oxff; 


if(abs((int)*input)<126)  /*  -1 28,-1 27,+1 27, 126  are  reserved  codes*/ 

{ 

*treeptr++  =  (signed  char)*input; 
outcount-i-i-; 

} 

else  if  (abs((int)*input)<=32895) 

{ 

*treeptr+-i-=-128; 

*treeptr+-i-=  (signed  char)((int)*input»8); 

*treeptr++=  (unsigned  char)((int)*input&lowbyte); 
outcount+=3; 

if(DEBUG)  printfC'decomp  %i  into  %i  +  %i\n",(int)*input, 
(int)*(treeptr-2),(int)*(treeptr-1)); 

} 

else 

{ 

printf("ABORTED:  range  error  exceeded  on  abs  calc  =  %i\n",(int)*input); 
fwrite(&treeptr[-outcount],sizeof(signed  char),outcount,outfile); 
printf("%i  token  saved  in  outfile\n'',outcount); 
exit(O); 

} 

return  treeptr; 


int  threshcheck(float  *locQ, int  level, int  depth, int  retval) 

{ 

int  i,  chkval; 

chkval  =  threshold1+(threshold2*depth); 

for  (i=0;i<depth;i++)  /‘left  side*/ 

{ 

if(abs(loc[level][i])>chkval)  retvall=2; 

} 

for  (i=depth;i<depth‘2;i++)  /*  right  side*/ 

{ 

if(abs(loc[ievel][i])>chkval)  retvail=1 ; 


} 

if  (retval==3)  return  retval;  /*  both  sides  signif,  stop  7 
else  if  (depth<lookahead&&levei>1 ) 

return  threshcheck(loc, level-1  ,depth+1 , retval): 
else  return  retval; 

} 


************★★*/ 


signed  char  *treetraverse(float  'vector, signed  char  'treeptr,  float  *loc[], 
float  *treeQ,int  level) 

{ 

int  i; 

int  nextievel; 
int  initval  =0; 
int  initdepth  =1 ; 


/'encode  left  node  and  subtree  if  not  zero'/ 


if  (abs('loc[level])<=threshold1&&level>1) 

{  /'  need  to  see  if  anything  sig  is  below'/ 

switch(threshcheck(ioc, level-1, initdepth, initval)) 

{ 

case  0:  /'  both  subtrees  insignificant '/ 
if  (DEBUG2) 

{ 

printf ("level  %i,  stop  at  %i,  i=%i,  o=126  -  both  pruned, " 
"going  back\n",  level, (int)(loc[level]-tree[level]), 

(int)  (signed  char)'loc[level]); 

} 

'treeptr-n-  =  126; 
counti  26-1-1-; 
outcount-n-; 

nextievel  =  level;  /'update  subtree  pointers  2  ea  level'/ 
for(i=nextlevel:i>0:i-) 

{ 

loc[i]  -H=  1«(nextlevel-i); 

} 

return  treeptr; 

case  1 :  /'  left  subtree  insignificant,  process  right'/ 
if  (DEBUG2) 

{ 

printf  ("level  %i,  right  at  %i,  i=%i,  o=127  -  pruned  left, " 
"going  right  \n",  level,(int)(loc[level]-tree[level]), 
(int)(signed  char)'loc[ievel]); 

} 

'treeptr-4-f  =  127; 
counti  27-H-: 
outcount-n-; 


/'update  subtree  pointers,  before  doing  right '/ 
loc[level]-i-i-;  /'  consume  input '/ 


nextievel  =  level-1: 
for(i=nextlevel:i>0:i-) 

{ 

loc[i]  +=  1«(nextlevel-i); 

} 

treeptr  =  treetraverse(vector, treeptr.loc, tree, level-1 ); 
return  treeptr; 

case  2:  /*  right  subtree  insignificant,  process  left*/ 
if  (DEBUG2) 

{ 

printf{“level  %i,  left  at  %i,  i=%i,  o=-127  -  pruned  right, " 
"going  left  \n",  level,(int)(loc[level]-tree[level]), 
(int)(signed  char)*loc[level]): 

} 

*treeptr-H+  =  -127; 

outcount-n-: 

countm127-i-H; 

treeptr  =  treetraverse(vector,treeptr,loc,tree,level-1): 

if  (DEBUG2)  printfC'level  %i,  right  at  %i  skipped\n", level, 
(int)(loc[level]-tree[level])); 

/"update  subtree  pointers,  after  left  done  */ 
nextievel  =  level-1 ; 
for(i=nextlevel;i>0:i-) 

{ 

loc[i]  -^=  1«(nextlevel-i); 

} 

loc[level]-i-i-: 
return  treeptr; 

/*  case  3:  both  sides  are  significant  fall  through  and  process  */ 

} 


/*  or  else  process  normally  */ 

treeptr  =  encode(treeptr,loc[level]); 

if(level>1) 

{ 


if  (DEBUG2)  printfC'level  %i,  left  at  %i,  i=%i,  o=%i\n", 
level,loc[level]-tree[leveli,(int)((signed  char)*loc[level]), 
(int)((signed  char)*(treeptr-1 ))); 

treeptr  =  treetraverse(vector,treeptr,loc,tree,level-1); 

/‘encode  right  subtree  if  level>1  */ 

/* 

if  (DEBUG2)  printfC'level  %i,  right  at  %i,  i=%i,  o=%i\n", 
level,locilevel]-tree[level],(int)((signed  char)*loc[level]), 
(int)((signed  char)*(treeptr-1 ))); 

*/ 


if  (DEBUG2)  printf("level  %i,  right  at  %i\n", level, 
loc[level]-tree[level]): 


treeptr  =  treetraverse(vector,treeptr,loc, tree, level-1): 

} 

else 


if  (DEBUG2)  printf("level  %i,  done  at  %i,  i=%i,  o=%i\n", 
level, loc[level]-tree[level],(int)((signed  char)*loc[level]), 
(int)((signed  char)*(treeptr-1 ))); 

} 

loc[level]-(-i-;  /*  update  location  */ 
return  treeptr; 

} 


void  zerotree(float  ‘vector, signed  char  *treevector,int  size,int  levels) 

{ 

int  i,  maxlevel; 
float  “loc; 
float  “tree; 
signed  char  ‘treeptr; 


printfC'setting  up  tree  pointers\n"); 
i=0; 

while(1«i<size)  maxlevel=+-i-i: 
loc  =  getvect(maxlevel+1,sizeof(int)); 
tree  =  getvect(maxievei+1,sizeof(int)): 

printfC'maxlevel  is  %i,  levels  requested  are  %i\n", maxlevel, levels); 

if(levels>maxlevel) 

{ 

printf("too  many  levels  requested\n''); 
exit(O); 

} 


tree[maxlevel]  =  vector;  !*  should  add  test  for  1>  2“levels  ‘/ 
loc[maxlevel]  =  vector; 

for  (i=0;i<maxlevel;i4-i-)  /‘  set  up  pointers  0  =  end,  maxlevel  =  start‘/ 

{ 

tree[i]  =  &vector[size»i] ; 
loc[l]  =  tree[i]; 

} 


if(DEBUG)printf("vector  address  =  %i\n",(int)vector); 

for  {i=maxlevel;i>0;i-) 

{ 

printfC  loc  %i  offset  %i  contains  %i  floats", 
i,loc[i]-loc[maxlevel],loc[i-1]-loc[i]): 


101 


if(DEBUG)  printf(",  starts  at  addr  %i'',(int)loc[i]); 
printf("\n"): 

} 

treeptr  =  treevector; 
i=0: 

while  (&vector[i]<  tree[levels]) 

{ 

treeptr  =  encode(treeptr,&vector[i]); 
if(DEBUG)  printf('Vect  =  %6.2f  gives  c  =%i  \n", 
vector[i],(int)treeptr[-1]): 
i++: 

} 

printf(“%i  tokens  copied\n",i): 

while  (loc[levels]<tree[levels-1])  /*  begin  recursive  tree  search  call  */ 

{ 

treeptr  =  treetraverse(vector, treeptr, loc, tree, levels); 

} 

} 


signed  char  *decode(signed  char  *treeptr,float  *outputloc) 
int  temp; 

if  (*treeptr==-128)  /*  rebuild  integer  value  */ 

{ 

treeptr++; 

temp  =  (int)((signed  char)*treeptr++); 
temp  «=8; 

temp  +=  (int)((unsigned  char)*treeptr++); 

*outputloc  =  (float)temp; 
outcount++; 

if  (DEBUG)  printfC'rebuild  %i  from  %i  +  %i\n", 

(int)*outputloc,(int)*(treeptr-2),(int)*(treeptr-1)); 

} 

else 

{ 

*outputloc=(float)*treeptr++; 

outcount++; 

} 

return  treeptr; 


signed  char  *buildtree(float  *tree[],  signed  char  *treevector,int  size, 
signed  char  ‘treeptr,  int  level,int  offset) 

{ 


signed  char  *treelast: 


if  (offset>=tree[level-1  ]-tree[level]) 

{ 

printf("level  %i  .size  %i,  index  %i  completedVn", 
level, tree[level-1]-tree[level], offset): 
return  treeptr; 

} 

switch  (*treeptr) 

{ 

case  126: 
if  (DEBUG2) 

{ 

printfC'level  %i,  stop  at  %i,  i=%i,  o=0  -  both  pruned,  going  back  \n", 
level,(int)offset,(int)(signed  char)*treeptr); 

} 

/*  consume  token  */ 
outcount++; 
return  ++treeptr; 


case  127: 
if  (DEBUG2) 

{ 

printfC'level  %i,  right  at  %i,  i=%i,  o=0  -  pruned  left,  going  right  \n", 
level, (int)offset,(int)(signed  char)*treeptr): 

} 

outcount++: 

treeptr++; 

/*if  input  exhasted  quit  7 
if  (offset>=tree[levei-1  ]-tree[level]) 

{ 

printfC'level  %i  .size  %i,  index  %i  completed\n", 
level,tree[level-1]-tree[level],offset); 
return  treeptr; 

} 


if(level>1 )  /*  continue  7 

{ 

treeptr  =  buildtree(tree,treevector,size, treeptr, 
level-1 ,2*(offset)+1); 

} 

return  treeptr; 


case -127:  /*skipright7 

if  (DEBUG2) 

{ 

printfC'level  %i,  left  at  %i,  i=%i,  o=0  -  prune  right,  going  left  \n", 
level, (int)offset,(int)(signed  char)*treeptr); 

} 

outcount++; 

treeptr-i-+; 


/*if  input  exhasted  quit  */ 
if  (offset>=tree[level-1]-tree[level]) 

{ 

printfC'level  %i  .size  %i,  index  %i  completed\n'', 
level, tree[level-1]-tree[level], offset); 
return  treeptr; 

} 

/*go  left  if  level>1  7 
if(level>1) 

{ 

treeptr  =  buildtree(tree,treevector, size, treeptr, level-1 , 
2*  (offset)); 

} 

if  (DEBUG2)  printfC'level  %i,  right  at  %i  skipped\n’', level, 
offset); 
return  treeptr; 


default:  /*  process  normally  7 

/*go  left  if  level>1 )  7 

treeptr  =  decode(treeptr,&tree[level][offset]); 

if(level>1) 

{ 


if  (DEBUG2)  printfC'level  %i,  left  at  %i,  i=%i,  o=%i\n", 
level, offset,(int)((signed  char)*(treeptr-1)), 
(int)((float)tree[level][offset])); 

treeptr  =  buildtree(tree,treevector,size,treeptr, 
level-1 ,2*(offset)); 

/*  and  then  do  the  right  side  7 

/*if  input  exhasted  quit  7 
if  (offset>=tree[level-1  ]-tree[level]) 

{ 

printfC'level  %i  .size  %i,  index  %i  completed\n", 
level, tree[level-1]-tree[level], offset); 
return  treeptr; 

} 

if  (DEBUG2)  printfC'level  %i,  right  at  %i\n", 
level, offset); 

treeptr  =  buildtree(tree,treevector,size,treeptr, 
level-1 ,2*(offset)-i-1); 

} 

else 

{ 

if  (DEBUG2)  printfC'level  %i,  done  at  %i,  i=%i,  o=%i\n". 


level, offset, (int)((signed  char)*(treeptr-1)), 
(int)((float)tree[level][offset])); 

} 

} 

return  treeptr; 

} 


void  unzerotree(float  *vector,signed  char  *treevector,int  size,int  levels) 

int  i,  maxlevel; 
float  **loc; 
float  **tree; 
signed  char  ‘treeptr; 
int  offset  =  0; 


printfC'rebuilding  tree  pointers\n"); 
i=0; 

while(1«i<size)  maxlevel=++i; 
loc  =  getvect(maxlevel+1  ,sizeof(int)); 
tree  =  getvect(maxlevel+1  ,sizeof(int)); 

printf("maxlevel  is  %i,  levels  to  restore  are  %i\n",maxlevel, levels); 

if(levels>=maxlevel) 

{ 

printf("too  many  levels  requested\n''); 
exit(O); 

} 


tree[maxlevel]  =  vector; 
locfmaxlevel]  =  vector; 

for  (i=0;i<maxlevel;i++)  /‘set  up  pointers  0  =  end,  level+1  =  begining  */  { 
tree[i]  =  &vector[size»i] ; 
loc[i]  =  tree[i]; 

} 


if  (DEBUG)  printf("vector  address  =  %i\n“,(int)vector); 
for  (i=maxlevel;i>0;i--) 

{ 

printf("  loc  %i  offset  %i  contains  %i  floats", 
i,loc[i]-loc[maxlevel],loc[i-1]-loc[i]); 
if  (DEBUG)  printfC,  starts  at  addr  %i'',(int)  loc[i]); 
printf("\n"); 

} 


treeptr  =  treevector; 
i=0; 

while  (&vector[i]<  tree[levels]) 


{ 

treeptr  =  decode(treeptr,&vector[i]): 
if(DEBUG)  printf("vect  =  %6.2f  \n",vector[i]); 
i++; 

} 

printf(”%i  tokens  copied\n'',i); 
offset=0; 

/*  while  there  is  still  input*/ 

while(offset<=tree[levels]-tree[maxlevel]) 

{ 

treeptr  =  buildtree(tree,treevector, size, treeptr, levels, offset++); 

} 


} 


void  main(int  argc,  char  *argv[]) 

{ 

int  fsize; 

int  xdim=0,  ydim=0,  size=0; 
float  ‘vector; 
float  ‘headvector; 
signed  char  ‘treevector; 
float  flevelq,  flevelzt; 

char  string[25]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,01; 

setoptions(argc,argv); 

infile  =  fopen(argv[1],"rb"); 
if(infile  ==  NULL) 

{ 

printf("File  open  failure\n"); 
exit(O); 

1 

strcpy(string,argv[1]); 


if  (direction>0)  strcat(string,".zt"); 
else  strcat(string,".tz“); 

outfile  =  fopen(string,"wb“); 
if(outfile==  NULL) 

{ 

printfC'File  open  failure\n''); 
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exit(O): 

} 

size  =  filesize(infile); 

if  (hsize==-1)  /*  no  header  specified  */ 

{ 

fread(&flevelq,sizeof(float),1, infile);  /‘levels  processed  by  quant*/ 
hsize  =  (flevelq+info)*2: 

} 

if(direction>0)  size  =  ((size-1  )/4)-hsize: 
else  size  -=  (hsize+info); 
head  vecto  r  =getvect(  hsize,  sizeof  (float) ) ; 
treevector  =getvect(size*2,sizeof(signed  char)); 
if  (direction==1) 

{ 

vector  =  getvect(size,sizeof(float)); 

fwrite(&flevelq,sizeof(float),info,outfile);/*levels  processed  by  quant*/ 
flevelzt  =  (float)level; 

fwrite(&flevelzt,sizeof(float),1,outfile);/*levels  process  byzt  */ 
fread(headvector,sizeof(float), hsize, infile);  /*get  data  vector*/ 
fwrite(&size,sizeof(float),1,outfile); 
printf(“saving  header  of  size  %i\n'', hsize); 
fwrite(headvector,sizeof(float),hsize,outfile);  /‘save  data  vector*/ 
fread(vector,sizeof(float),size,infile);  /‘get  image  vector  */ 
zerotree(vector,treevector,size,level); 
fwrite(treevector,sizeof  (signed  char),outcount,outfile); 
printf("%s  =  %i  tokens,  converted  to  %i  char\n",argv[1],size,outcount); 
printfC'output  126  =  %i,  127  =  %i,  -127  =  %i  total  code  tokens  %i\n", 
count126,count127,countm127,count126-t-count127+countm127); 
printf("token  overhead  %3.2f\n", 

1 00.0*((float)(count1 26-t-countl  27+countm127)/(float)outcount)); 

} 

else 

{ 

fread(&flevelzt,sizeof(float),1,infile);/*level  processed  by  zt  */ 

level  =  (int)flevelzt; 

fread(&fsize,sizeof(float),1, infile); 

vector  =  getvect(fsize,sizeof(float)); 

fwrite(&flevelq,sizeof(float),info,outfile); 

fread(headvector,sizeof(float),hsize,infile);  /‘get  data  vector*/ 

printf("saving  header  of  size  %i\n",hsize); 

Write(headvector,sizeof(float),hsize,outfile);  /‘save  data  vector*/ 
fread(treevector,sizeof(signed  char),size,infile); 
unzerotree(vector, treevector, fsize, level); 
fwrite(vector,sizeof(float),fsize,outfile); 

printf("%s  =  %i  tokens,  convert  to  %i  floats,  expanded  to  %i\n",argv[1],size,outcount,fsize); 

} 

free(vector); 

free(treevector); 

fclose(infile); 

fclose(outfile); 

} 


/*  - . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  WAVE.C 

PURPOSE:  Perform  a  linear  orthogonal  transform  on  a  floating  point  data  file 

Input  file  is  a  binary  float  file 
output  file  is  a  float  file 

FEATURES:  Has  several  options  to  control  the  selection  of  filters  and  the 
number  of  levels  processed,  allows  snaking  of  the  input  on  the  row  length 
to  minimize  the  edge  effect  of  the  image  on  the  high  band  pass 


•*/ 


#include<stdio.h> 

#include<stdlib.h> 

#include<string.h> 

#define  MAXSIZE  262145 


GLOBALS 


*******************★**********:*************************★*** 


struct  controlblock  { 
char  infile[81]; 
char  outfile[81]; 
char  root[81]; 
char  extension[5]; 
int  filter; 
int  compress; 
char  action; 
int  level; 
int  shift; 
double  normal; 
int  rotate; 
int  snake; 

}: 

/*  daubechies  filters  2  -  20  7 
double  filtcoef[1 1][20]  = 

{ 

{ 1.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/*  daub  2  7 

{  0.707106781187,  0.707106781 187,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 


0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/*  daub  4  */ 

{ 0.482962913145,  0.836516303738,  0.224143868042,-0.129409522551, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/*  daub  6  */ 

{ 0.332670552950,  0.806891509311, 0.459877502118,-0.135011020010, 
-0.085441273882,  0.035226291882,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/*  daub  8  */ 

{ 0.230377813309,  0.714846570553,  0.630880767930,-0.027983769417, 
-0.187034811719,  0.030841381836,  0.032883011667,-0.010597401785, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/*  daub  10  V 

{  0.160102397974,  0.603829269797,  0.724308528438,  0.138428145901, 
-0.242294887066,-0.032244869585,  0.077571 493840,-0.006241 49021 3, 
-0.012580751999,  0.003335725285,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/‘daub  12  */ 

{  0.111540743350,  0.494623890398,  0.751133908021,  0.315250351709, 
-0.226264693965,-0. 1 29766867567,  0.097501 605587,  0.027522865530, 
-0.031 58203931 8,  0.000553842201 , 0.00477725751 1  ,-0.001 077301 085, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/*  daub  14*/ 

{  0.077852054085,  0.396539319482,  0.729132090846,  0.469782287405, 
-0.143906003929,-0.224036184994,  0.071309219267,  0.080612609151, 
-0.038029936935,-0.01 6574541 631 , 0.01 2550998556,  0.000429577973, 
-0.001801640704,  0.000353713800,  0.000000000000,  0.000000000000, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/‘daub  16*/ 

{  0.054415842243,  0.312871590914,  0.675630736297,  0.585354683654, 
-0.01 58291 05256,-0.28401 5542962,  0.000472484574,  0.1 28747426620, 
-0.01 7369301X)02,-0.044088253931 ,  0.01 3981 02791 7,  0.008746094047, 
-0.004870352993,-0.000391 740373,  0.000675449406,-0.0001 1 7476784, 
0.000000000000,  0.000000000000,  0.000000000000,  0.000000000000  }, 
/*  daub  18  */ 

{  0.038077947364,  0.243834674613,  0.604823123690,  0.657288078051, 
0.133197385825,-0.293273783279,-0.096840783223,  0.148540749338, 
0.030725681479,-0.067632829061, 0.000250947115,  0.022361662124, 
-0.004723204758,-0.004281 503682,  0.001 847646883,  0.000230385764, 
-0.000251963189,  0.000039347320,  0.000000000000,  0.000000000000  }, 
/*  daub  20  */ 

{  0.026670057901, 0.188176800078,  0.527201188932,  0.688459039454, 
0.281 172343661  ,-0.249846424327,-0.195946274377,  0.127369340336, 


0.093057364604,-0.071 3941 471 66,-0.029457536822,  0.03321 2674059, 
0.003606553567,-0.010733175483,  0.001395351747,  0.001992405295, 
-0.000685856695,-0.0001 16466855,  0.000093588670,-0.000013264203  } 

}: 


/*  openfile  -  open  file  and  check  for  errors  */ 


FILE  *openfile{char  *file,  char  *opt) 

{ 

FILE  *ptr: 

if((ptr  =  fopen(file,  opt))  ==  NULL) 

{ 

fprintf(stderr,  "Can  not  open  the  file  %s\n",  file); 
exit(1); 

} 

return(ptr); 

} 


***************************  / 


/*  loadfile  -  loads  in  a  data  file  of  floating  point  numbers  7 


★****★★*★★****/ 


long  loadfile(float  ‘data,  struct  controlblock  control) 

{ 

FILE  *in; 
int  i,j; 
long  size; 

int  width,dblwidth,hlfwidth,span; 
long  len; 
float  tmp; 

/*  open  the  input  file  */ 
printf("loading  %s\n'',control.infile): 

{ 

in  =  openfile(control.infile,  "rb"); 

fseek(in,  OL,  2); 

len  =  ftell(in); 

fseek(in,  OL,  0); 

if(len%sizeof(float)) 

{ 

fprintf(stderr,  "Error  in  reading  input  file.\n"); 
exit(1); 

} 

else 

{ 

size  =  len/sizeof(float): 
if(size  >  MAXSIZE) 

{ 


fprintf(stderr,  "Input  file  is  too  largeAn"); 
exit(1): 

} 

else 

{ 

if(fread(data,  sizeof(float),  size,  in)  !=  size) 

{ 

fprintf(stderr,  "Error  reading  input  fileAn"); 
exit(1): 

} 

} 

} 

if  (control.action!='i‘&&control.snake>0)  /*  reverse  every  other  row  V 

{ 

width  =  control.snake; 

printfC'snake  option  with  width  %i  used\n", width); 
dbiwidth  =  width*2; 
hlfwidth  =  width/2; 
for  G=width;j<size;j+=dblwidth) 

{ 

span  =  j+width-1 ; 
for(i=0;i<hlfwidth;i++) 

{ 

tmp  =  dataO+i]; 
dataO+i]  =  data[span-i]; 
data[span-i]  =  tmp; 

} 

} 

} 

} 

fclose(in); 
return  (size); 

} 


/*  extendedhex  -  This  routine  takes  an  integer  between  0  and  36  and  converts 
it  into  a  character  that  goes  from  0-9, A-Z.  */ 


char  extendedhex(int  value) 

{ 

if{(value  >  35)  II  (value  <  0)) 

{ 

fprintf(stderr,  “Value  too  big  to  specify  the  file  extension  properlyAn"); 
exit(1); 

} 

if(value  <  10)  return(’0'  +  value); 
else  return('a'  +  value  - 10); 

} 


r  savefile  -  This  routine  saves  results  in  output  using  the  coined  extension  */ 


^***-k**ifk**-kicie-kie**-kirie'k1fkieie*ickifif**iritie1c***ie**ic**1ci(iti(i(ick-k'k*-k*  j 


void  savefile(float  *data,int  level, long  length, struct  controlblock  control) 

{ 


FILE  *out: 
int  i,j,size; 

int  width,dblwidth,hlfwidth,span; 
float  tmp; 


/*  construct  the  output  name  */ 
if(strlen(control.root)  !=  NULL) 

{ 

control.extension[0]  = 

r  use  the  action  type  for  the  first  byte  7 
if(control. action  ==  's') 

{ 

control. extension[1]  =  extendedhex(control.level); 

} 

else 

{ 

control. extension[1]  =  control.action; 

} 


/*  use  the  wavelet  type  for  the  second  byte  7 
control.extension[2]  =  extendedhex(2*control.filter); 

/*  use  the  level  of  the  decomposition  for  the  third  byte  7 
control.extension[3]  =  extendedhex(level): 
control.extension[4]  =  0x00; 
strcpy(control.outfile,  control. root); 
strcat(control.outfile,  controi.extension); 

} 

if  (control.action=='i'&&control.snake>0)  /*  reverse  every  other  row  7 

{ 

width  =  control. snake; 

printf ("snake  option  with  width  %i  usedNn”, width); 
dbiwidth  =  width*2; 
hlfwidth  =  width/2; 
for  (j=width;j<length;j+=dblwidth) 

{ 

span  =  j+width-1; 
for(i=0;i<hlfwidth;i++) 

{ 

tmp  =  dataO+i]; 
data[j+i]  =  data[span-i]; 
data[span-i]  =  tmp; 

} 

} 

} 


/*  open  the  output  file  and  save  7 
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printf("saving  %s\n", control. outfile); 
out  =  openfile(control. outfile,  "wb"); 
if{fwrite(data,  sizeof(float),  length,  out)  !=  length) 
{ 

fprintf(stderr,  "Error  writing  ouput  file.\n"); 
exit(1): 

} 

fclose(out); 

} 


^**************★**************★*★****★**********★******1^**  j 

r  getvector  -  memory  allocation,  error  check  and  initialize  7 


void  *getvector(int  size,  long  length) 

{ 

char  *ptr; 
long  i; 

ptr  =  (char  *)mailoc(size*length): 
if(ptr  ==  NULL) 

{ 

fprintf(stderr,  "Error  in  allocating  memory.\n“): 
exit(1): 

} 

/*  initialize  7 

for(i=0;  i<size*length:  i++)  ptr[i]  =  0; 
return((void  *)ptr); 

} 


/*  copyvect  -  copy  one  vector  to  another,  pad  left  on  positive,  right 
on  negative  7 


void  copyvect(float  *out,  float  *in,  long  length,  int  overrun) 

{ 

long  i; 

printfC'vector  copy,  size  %i\n“,length); 

if(overrun  <  0) 

{ 

overrun  *=  -1 ; 

for(i=0;  i<length;  i++)  out{overrun  +  i]  =  in[i]; 

r  i%length  allows  multiple  repetitions  if  input  shorter  than  overrun  7 
for(i=0:  i<overrun;  i++) 

out[overrun-1 -i]  =  out[overrun+length-1-i%length]: 

} 

else 

{ 


for(i=0:  i<length;  i++)  out[i]  =  in[i]: 

r  continually  repeat  while  necessary  V 

for(i=0:  i<overrun;  i++)  out[length+i]  =  out[i%length]; 

} 

} 


/*  minsize  -  calculate  a  minimum  length  that  is  power  of  2  and 
wil  contain  the  data  V 


long  minsize(int  ‘iterations,  long  number) 

{ 

int  loops  =  0; 
long  pow2  =  1 ; 

while(1) 

{ 

if(pow2  >=  number) 
break; 
pow2  «=  1 ; 
loops++; 

} 

r  the  analysis  can’t  exceed  the  length  of  the  data  */ 
if(*iterations  ==  0)  ‘iterations  =  loops; 
else  if(‘iterations  >  loops) 

{ 

fprintf(stderr,  "Data  not  long  enough  for  the  number  of  iterations  requested\n"); 

fprintf(stderr,  "Max  #  of  iterations  possible  is  %d.\n“,  loops); 

exit(1); 

} 

return  (pow2); 

} 


/*  buildfilter  -  sets  up  wavelet  coefficients  for  transform, 
if  filt  is  negative  then  sets  up  reconstruction  coefficients  ‘/ 


void  buildfilter(dbuble  ‘even,  double  ‘odd,  int  filter) 
int  i; 

if(filter  >  0) 

{ 

for{i=0;  i<2‘filter;  i++) 

{ 

r  low  pass  ‘/ 
even[i]  =  filtcoef[filter][i]; 

/‘  high  pass  ‘/ 

odd[i]  =  filtcoef[filter][2*filter  -  i  -1]; 
if(i%2)  odd[i]‘=-1; 
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} 

} 

else 

{ 

filter  *=-1; 
for(i=0:  i<filter:  i++) 

{ 

/*  inverse  even  index  */ 

even[i*2]  =  filtcoef[filter][filter*2  -  2*(i+1)]; 

even[i*2+1]  =  filtcoef[filter][2*i+1]; 

/*  inverse  odd  index  */ 

odd[i*2]  =  filtcoef[filter][filter*2  -  (2*i+1)]; 

odd[i*2+1]  =  -filtcoef[filter][i*2]; 

} 

} 

} 


/*  wavelet  -  performs  wavelet  decomposition,  returns  high  frequency  component 
in  first  half  of  the  vector  and  low  frequency  in  last  half  of  the  vector.  */ 


void  wavelet(float  *buf,  long  length,  int  width,  double  *evenfilter, 
double  *oddfilter) 

{ 

long  i; 
intj: 

int  times; 
float  even,  odd; 


printfC'transforming,  width  %i\n",width); 

/*  begin  the  loop  through  the  data  */ 
for(i=0;  i<length;  i  +=  2) 

{ 

even  =  0.; 
odd  =  0.; 

for{j=0;  j<width;  j++) 

{ 

even  +=  evenfilter|j]*buf[i+j]; 
odd  +=  oddfilter|j]*buf[i+j]; 

} 

buf[i]  =  even; 
buf[i+1]  =  odd; 

} 

} 


^★★★★★★★★★★★★★★★★^■fci^t**************************^^***********  j 

r  decimate  -  permutes  the  vector  into  a  second,  decimate  if  length  positive 


or  interpolate  if  negative,  neg  overrun  right  to  left,  pos  overrun  left  to 
right  7 


void  decimate(float  *out,  float  *in,  long  length,  int  overrun) 

{ 

long  i; 

long  halflength; 
int  offset; 


printfC'vector  permute,  length  %i\n“,length); 

if(length  <  0)  halflength  =  -length/2; 
else  halflength  =  length/2; 

if(overrun  <  0)  offset  =  -overrun; 
else  offset  =  0; 

if(length  <  0) 

{ 

/*  decimate  7 
length  *=  -1 ; 

for(i=0:  i<halflength;  i+-^) 

{ 

r  copy  the  even  index  -  low  pass  result  7 
out[offset+i*2]  =  in[i]; 

/*  copy  the  odd  index  -  high  pass  result  7 
out[offset-i-i*2  -I- 1]  =  in[i+half length]; 

} 

} 

else 

{ 

r  split  the  alternate  values  up  7 
for(i=0;  i<halflength;  i-n-) 

{ 

outfoffset+i]  =  in[i*2]; 
out[offset+i-Hhalflength]  =  in[i*2+1]; 

} 

} 

if(overrun  <  0)  - 

{ 

/*  pad  the  left  side  7 
overrun  *=  -1; 

for(i=0;i<overrun;i++) 

out[overrun-1-i]  =  out[overrun-i-length-1-i%length]; 

else 

{ 

/*  pad  the  right  side  7 

for(i=0;  i<overrun;  i-n-)  out[length+i]  =  out(i%length]: 
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} 


} 


y***************************************-*r  ***************** ! 

r  rotate  -  shifts  the  vector  right  or  left  based  on  the  value  of  dir, 
left-over  data  values  are  wrapped  to  the  other  end  7 
^★★★★★★★★★★★★★★**********^****************************^***^ 


void  rotate(float  *data,  long  length,  int  shift,  int  dir) 

{ 

float  bfs[20];  /*  is  good  up  to  D40  filter  7 
long  i; 

r  if  shift  distance  exceeds  the  length  then  reduce  the  shift  7 
shift  %=  length; 

if(dir  ==  1 )  /*  rotate  upward  7 

{ 

for(i=0;i<shift;i++)  bfs[i]  =  data[length-shift+i]: 

for(i=0:i<length-shift:i++)  data[length-i-1]  =  data[length-shift-i-1]; 
for(i=0;i<shift;i-i-i-)  data[i]  =  bfs[i]: 

} 

if(dir  ==  -1 )  /*  rotate  downward  7 

{ 

for(i=0;i<shift:i+-i-)  bfs[i]  =  datafi]; 

for(i=0;i<length-shift;i++)  data[i]  =  data[i+shift]; 
for(i=0;i<shift:i++)  data[length-shift+i]  =  bfs[i]; 

} 

} 


^★★********************************************jif******ik*** ! 

r  help  -  prints  options  list  to  the  screen  7 

void  help(void) 

{ 

printf("wavelet  <input>  (options)\t  default  options:  -d4  -n2  -t\n"); 
printf("Options:\n"); 

printf(“-o  <output>\n\tOutput  file  root  or  whole  nameAn"); 

printf(“-d#\tSpecify  the  wavelet  transform  d2  through  d20An"): 

printf(''-n#\t1/2  Normalization.  1-On  inverse,  2-Root  on  both,  3-On  transform.\n"); 

printf(''-t#\tPerform  the  transform  to  #  levels  .\n''); 

printf(“-i#\tPerform  the  inverse  transform  to  #  levels  of  the  pyramid.\n“); 

printf("-w#\tSnake  input  data  on  row  width  given.  \n"); 

printf(“-s#\tSame  as  -t#  but  save  each  level  as  a  separate  flle.Xn"); 

printf("-l#\tPerform  a  low  pass  filter  on  the  data  transformed  to  the  #  level.\n"); 

printf(''-h#\tPerform  a  high  pass  filter  on  the  data  transformed  to  the  #  level.Nn"): 

printf(''\t(Analyze  to  the  highest  level  possible  if  #  is  not  specified.)\n"); 

printf("-r\tCorrelate  coefficients  at  different  levels  by  shifting  results.\n“): 

} 

^******************************^*****************^********^ 

/*  setcontrol  -  sets  command  line  options  in  'controiblock'7 


^iticifk-kic-k*-kit**it**it***’k-k**ie-kitifk-k**-kieifie-k*it*ifk**-k*itifk1tie-kitieie'k**  J 

Struct  controlblock  setcontrol(int  argc,  char  *argv[]) 

{ 

int  i; 

int  level; 
char  *ptr; 

struct  controlblock  control; 

/*  set  up  the  default  values  7 
strcpy(control.infile,  "wavelet.dat"); 
strcpy(control.outfile, 
strcpy(control.root, 
control.filter  =2; 
control.compress  =  0; 
control. normal  =  1 .0; 
control.action  =  't'; 
control. level  =  0; 
control.shift  =  0; 
control,  rotate  =  0; 
control.snake  =  0; 

for(i=1 ;  i<argc;  i++) 

{ 

if(argv[i][0]  == 

{ 

switch(argv[i]il]) 

{ 

case 

helpO; 
exit(O); 
break; 
case  'o': 

strcpy(control .  outf  i  le,  arg v[i+1  ]) ; 
i++; 
break; 
case  '1': 

control.shift  =  1 ; 
break; 
case  'n': 

level  =  atoi(&argv[i][2]); 
if  ((level<1)ll(level>3)) 

{ 

fprintf(stderr,  “%s  is  not  a  valid  option.\n",  argv[i]); 
exit(1); 

} 

if  (atoi(&argv[i][2])!=2) 

control.normal  =  1.41421356237309504880; 
if  (atoi(&argv[i][2])>2)  control.normal  =  1/control. normal; 
break; 
case  'r': 

control. rotate  =  1; 
break; 
case 'd'; 


control.filter  =  atoi(&argv[i][2]); 
if((control.filter%2  ==  1)  II  (control.filter  <  2)  II 
(control.filter  >  20)) 

{ 

fprintf(stderr,  "%s  is  not  a  valid  option.\n",  argv[i]); 
exit(1); 

} 

else  control.filter  /=  2; 
break; 
case  'c'; 

control.compress  =  1; 
break; 
case  T; 

control.action  =  'I'; 
level  =  atoi(&argv[i][2]); 
if  (level>0)  control. level  =  level; 
break; 
case  'h'; 

control.action  =  'h'; 
level  =  atoi(&argv[i][2]); 
if  (level>0)  control. level  =  level; 
break; 
case  'b': 

control.action  =  'b'; 
level  =  atoi(&argv[i][2]); 
if  (level>0)  control. level  =  level; 
break; 
case  'k': 

control.action  =  'k'; 
level  =  atoi(&argv[i][2]); 
if  (ievel>0)  control. level  =  level; 
break; 
case  Y: 

control.action  =  T; 
level  =  atoi(&argv[i][2]); 
if  (level>0)  control. level  =  level; 
break; 
case 't': 

control.action  =  't'; 
level  =  atoi(&argv[i][2]); 
if  (level>0&&level<22)  control. level  =  level; 
break; 
case  's': 

control.action  =  's'; 
level  =  atoi(&argv[i][2]); 
if  (level>0)  control. level  =  level; 
break; 
case  'w': 

control.snake  =  atoi(&argv[i][2]): 
break: 

default :  fprintf(stderr,  "invalid  option  =  %d\n“,  argv[i]); 
exit(1); 
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else 

{ 

strcpy(control.infile,  argv[i]); 

} 

} 

/*  set  up  the  root  and  extension  of  the  output  files  7 
if(strlen(control.outfile)  ==  0) 

{ 

strcpy(control.outfile,  control. infile); 
ptr  =  strchr(control.outfile, 
if(ptr) 

ptrfO]  =  0x00; 

strcpy(control.root,  control. outfile); 

} 

else 

{ 

if(strchr(control. outfile, '.')  ==  NULL)  strcpy(control.root,  control.outfile); 
else  if(control.action  ==  ‘s') 

{ 

fprintf(stderr,  "No  output  file  extension  allowed  with  option  -s\n"); 


return(control); 

} 


int  main(int  argc,  char  *argv[]) 

{ 

long  i; 
intj; 

float  ‘data; 
float  *buf; 
long  num; 
long  pow2; 
long  width; 
long  offset; 

struct  controlblock  control; 
double  evenfilfer[20]; 
double  oddfilter[20]; 

int  wvshft[1 1][2]  =  {{0,0},{0,0},{0,1  },{0,2},{0,3},{1 ,3}, 
{1,4}.{1,5},{1,6},{1.7},{1.8}}; 

/*  get  command  line  options  */ 
if(argc  <  2) 

{ 

fprint("\nwavelet  <input  file>  -(options)\n"); 
fprint( "  Use  -?  for  a  help  file.\n''); 
exit(1); 

} 
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else 

{ 

control  =  setcontrol(argc,  argv); 

} 

/*  allocate  the  memory  V 

printf ("starting  program\n"); 

data  =  (float  *)getvector(sizeof (float),  MAXSIZE); 

buf  =  (float  *)getvector(sizeof(float),  MAXSIZE  +  2*(control.filter-1)); 

r  load  data  V 

num  =  loadfile(data,  control); 

r  determine  the  largest  power  of  two  that  covers  the  data  */ 
pow2  =  minsize(&control.level,  num); 

printfC'depth  =  %i,  for  file  size  %i,  %i  levels\n",pow2,num,control.level); 
/*  set  up  the  filter  */ 

buildfilter(evenfilter,  oddfilter,  control.filter); 

/*  apply  normalization  */ 
for(j=0;  j<control.filter*2;  j++) 

{ 

evenfilter[j]  *=  control. normal; 
oddfilterQ]  *=  control. normal; 

} 

/*  if  -1  is  used  shift  the  data  one  to  the  left  7 

if((control.shift==1)&&(control.action!='i'))  rotate(data,  pow2,  1,-1); 

/*  perform  filter  pass,  save  the  high  frequency  component  in  lower 
half  of  vector  and  the  low  frequency  component  in  upper  half  7 

width  =  pow2; 

for(j=0;  j<control. level;  j++) 

{ 

/*  don't  decompose  if  only  doing  the  inverse  transform  7 
if(control.action  !=  'i') 

{ 

/*  do  the  transform  7 

copyvect(buf,  data,  width,  2*(control.filter-1)); 
wavelet(buf,  width,  2*control.filter,  evenfilter,  oddfilter); 
decimate(data,  buf,  width,  0); 

/*  rotate  the  low  and  high  frequency  numbers  to  register  properly  7 
if  (control,  rotate) 

{ 

rotate(data,  width/2,  wvshft[control.filter][0],  1); 
rotate(data+width/2,  width/2,  wvshft[control.filter][1],  1); 

} 

} 


/*  shift  to  the  upper  half  so  the  decomposition  can  be  done  again  7 
width  /=  2; 

} 

/*  if  only  transform  requested  save  and  quit  7 
if(control.action  ==  't') 


{ 

savefile(data,  control.level,  pow2,  control); 
exit(O); 

} 

if(control.action  ==  ’s') 

{ 

savefile(data,  0,  width, control); 

offset  =  width; 

forG=0;  j<control.level;  j++) 

{ 

savefile(data  +  offset,  j+1,  offset,  control); 
offset  *=  2; 

} 

exit(O); 

} 

/*  if  filtering  remove  either  the  high  or  low  frequency  component  */ 
if((control. action  ==  'h')  II  (control. action  ==  'b')) 

{ 

for(i=0;  i<width;  i++)  data[i]  =  0.; 

} 

if(control.action  ==  'b') 

{ 

for(i=width*2;  i<pow2;  i++)  data[i]  =  0.; 

} 

if(control.action  ==  'k') 

{ 

for(i=width;  i<width*2;  i++) 
data[i]  =  0.; 

} 

if(control.action  ==  T) 

{ 

for(i=width;  i<pow2;  i++) 
data[i]  =  0.; 

} 

/*  set  up  the  reconstruction  filter  */ 
buildfilter(evenfilter,  oddfilter,  -control.filter); 

/*  apply  normalization  */ 
for(j=0;  j<control.filter*2;  j++) 

{ 

evenfilterO]  /=  control. normal; 
oddfilterjj]  /=  control. normal; 

} 

/*  reconstruct  the  input  from  the  transform  results  7 
for(j=0;  j<control. level;  j++) 

{ 
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/*  expand  to  the  next  largest  scale  to  decompose  again  */ 
width  *=  2; 

/*  shift  the  low  and  high  frequency  vectors  for  reconstruction  V 

if(control. rotate) 

{ 

rotate(data,  width/2,  wvshft[control.filter][0],  -1); 
rotate(data+width/2,  width/2,  wvshft[control.filter][1],  -1); 

} 


/*  perform  the  inverse  transform  7 
decimate(buf,  data,  -width,  -2*(control.filter-1)): 
wavelet(buf,  width,  2*control.filter,  evenfilter,  oddfilter); 
copyvect(data,  buf,  width,  0); 

} 


/*  if  -1  is  used  shift  the  data  one  to  the  left  7 
iffcontrol.shift  ==  1)  rotate(data,  pow2, 1,1); 

r  output  result  7 

savefile(data,  control. level,  num,  control); 
printf("program  complete\n"); 
return  0; 

} 


/* . . . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  QONLY.C 

PURPOSE;  Quantizes  an  input  float  vector  from  a  wavelet 
transform  file 

Input  file  is  a  binary  float  file 
output  file  is  a  char  file 

Creates  a  data  header  of  1  +  2*levels  processed  (floats) 

Containing  levels  processed  and  mean  and  standard  deviation  of  each 
level  processed. 

FEATURES;  Uses  a  bin  calculation  routine  that  is  based  on  3  times  the 
standard  deviation,  puts  93%  of  the  coefficients  into  64  bins  (+/-32) 


#include  <stdio.h> 

#include  <stdlib.h> 

#include  <math.h> 

#include  <fcntl.h> 

#include  <string.h> 

#define  DEBUG  1 

FILE  *infile,  *outfile,  *binfile; 


void  *  getvect(int  n,int  size) 

{ 

void  *vectptr; 


vectptr  =  (void  *)  calloc(n,size); 
if(vectptr  ==  NULL)  { 
printf("Error:  dynamic  memory  failVn"); 
exit(O); 

} 

return(vectptr): 

} 


^★♦★★♦★★★★★★★★★★★★★★★★imr****************************^ 

void  quant(float  *vector,float  *headvector,int  size, 
int  stddev,int  mean,int  debug) 

{ 

float  fstddev,  fmean,  frtsdv; 
int  i=0; 


fstddev  =  (float)stddev/100; 

frtsdv  =  (fstddev>1)  ?  (float)sqrt((double)fstddev):1; 
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fmean  =  (float)  mean; 
mean  *=100; 

headvector[0]  =  mean; 
headvector[1]  =  stddev; 

while  (i<size) 

{ 

vector[i]  =  (vector[i]-fmean)/frtsdv; 
i++; 

} 

} 


void  unquant{float  ‘vector, float  *headvector,int  size.int  debug) 

int  1=0,  stddev,  mean=0,  rebuilt; 
float  fstddev,  fmean,  frtsdv,  temp; 


mean  =  headvector[0]; 
stddev  =  headvector[1]; 


fstddev  =  (float)  (stddev)/1 00.0; 

frtsdv  =  (fstddev>1)  ?  (float)sqrt((double)fstddev):1; 

fmean  =  (float)  mean/100; 

printfC'mean  =  %6.2f,  stddev  =  %6.2f,  root  =  %6.2f\n",fmean,fstddev,frtsdv); 

while  (i<size) 

{ 

vector[i]=(vector{i]*frtsdv)+fmean; 

i++; 

} 

} 


void  sectquant(float  *vector,float  ‘headvector,  int  size,int  zeroout,int  quantize,int  debug) 
int  i; 

double  high=-1 0000000.0,  low=1 0000000.0; 
int  zeros=0,  zrun=0,  zcount=0,  outsiders=0; 

double  sum=0.0,  squares=0.0,  mean=0.0,  var=0.0,  stddev=0.0,  sd4=0.0; 
double  rtsdv=0.0; 

if  (quantize==-1)  unquant(vector,headvector,size,debug); 

for  (i=0;i<size;i++) 

{ 

high  =  (vector[i]>high)  ?  vector[i] :  high; 
low  =  (vector[i]<low )  ?  vector[i] :  low  ; 
sum  +=  vector[i]; 

squares  +=  (double)vector[i]*vector[i]/(double)size; 


} 

mean  =  sum/(double)size: 
var  =  squares  -  (mean*mean): 
stddev  =  sqrt(var); 
sd4  =  3*stddev; 

rtsdv  =  sqrt(stddev);  /*  zero  threshold,  sqrt  of  std  dev  */ 


for  (i=0;i<size;i++) 

{ 

if  ((double)abs(vector[i])>sd4)  outsiders++: 
if  ((double)abs(vector[i])<rtsdv) 

{ 

zeros++; 

if  (quantllzeroout)  vector[i]=0.0: 

} 

else  if(/*(double)abs(vector[i])<sd4&&7zeroout) 

{ 

zeros++; 

vector[i]=0.0; 

} 

if  ((double)abs(vector[i])<rtsdv) 

{ 

zcount++; 

} 

else  zcount  =  0; 

zrun  =  (zcount>zrun)  ?  zcount :  zrun; 

} 

printfC  high=  %6.2f,  low=  %6.2f,  slze=  %8i  \n“,  high,  low, size); 
printfC  mean  =  %4.1f,  var  =  %8.2f.  stddev  =  %6.2f\n",  mean,  var,  stddev); 
printf("  zerothresh  =  %6.2f,  zeros  =  %i,  longest  zero  run  =  %i\n", 
rtsdv,zeros,zrun); 

printfC  outsiderthresh  =  %6.2f,  outsiders  =  %i\n",sd4,outsiders): 


} 


if  (quantize==1 )  quant(vector,headvector,size,(int)stddev*100, 
(int)mean,debug); 


void  main(int  argc,  char  *argvQ) 

{ 

int  i,  quantize,  quant,  zero; 

int  xdim=0,  ydim=0,  size=0,  level=0,  sect=0,  cursize=0; 
float  ‘vector; 
float  ‘headvector; 
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char  string[25]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},  *suf=''.bq"; 


if(argc  <  5) 

{ 

printf ("Usage  pgm  <filename>  xdim  ydim  Q  Z  {1/-1}  \n"); 
prIntfCFirst  name  is  source  file\n"); 


printf("Q=Levels  to  quantize  \n"); 

printf("Z=Levels  to  zero  \n"): 
printf("1  quantize,  -1  unquantize,  0  none\n"); 
exit(O): 

} 

infile  =  fopen(argv[1],"rb''); 
if(infile  ==  NULL) 

{ 

printfC'File  open  failure\n"); 
exit(O); 

} 

strcpy(string,argv[1]): 
strcat(string,".q"); 
outfile  =  fopen(string,"wb"); 
if(infile  ==  NULL) 

{ 

printfC'File  open  failureVn"); 
exit(O); 

} 

xdim  =  atol(argv[2]): 
ydim  =  atoi(argv[3]); 
quant  =  atol(argv[4]): 
zero  =  atoi(argv[5]); 
quantize  =  atoi(argv[6]): 
level  =  quant+zero; 

if  (abs(quantize)>1)  quantize  =  1; 
size  =  xdim*ydim; 


vector  =  getvect(size,sizeof(float)): 
headvector=getvect(2*level,sizeof(float)); 

if  (quantize==1)  fread(vector,sizeof{float),size,infile): 

else 

{ 

fread{&level,sizeof(int),1, infile); 
fread(headvector,sizeof(float),2*level, infile); 
fread(vector,sizeof(float),size,infile); 

} 

for  (l=level;i>0;i--) 

{ 

cursize  =  size/(int)pow(2.0,(double)i); 
if  (level==i) 

{ 

printf("For  level  %i,  sect  %i  \n",  level,  1); 
sectquant(vector,headvector,cursize,0,quantize,1); 

} 

if(i<=zero)  for  (sect=2;sect<=2;sect++) 

{ 

prlntf("For  level  %i,  sect  %i  zeroing\n",i,sect); 
sectquant(&vector[cursize],&headvector[i*2], 
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cursize,1, quantize,!); 

} 

else  for  (sect=2:sect<=2;sect++) 

{ 

printfC'For  level  %i,  sect  %i  quantizing\n",i,sect); 
sectquant(&vector[cursize],&headvector[i*2], 
cursize,0,quantize,1): 

} 

} 


if  (quantize==1) 

{ 

fwrite(&level,sizeof(int),1,outfile); 

fwrite(headvector,sizeof(float),lever2,outfile); 

fwrite(vector,sizeof(float),size,outfile); 

free(vector); 

free(headvector); 

} 

else 

{ 

free(headvector); 

fwrite(vector,sizeof(float),size,outfile); 

free{vector); 

} 

fclose(infile): 

fclose(outfile); 

} 
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/*- . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME;  QONLY2.C 

PURPOSE:  Quantizes  an  input  float  vector  from  a  wavelet 
transform  file 

Input  file  is  a  binary  float  file 
output  file  is  a  char  file 

Creates  a  data  header  of  1  +  2*levels  processed  (floats) 

Containing  levels  processed  and  mean  and  standard  deviation  of  each 
level  processed. 

FEATURES:  Uses  a  bin  calculation  routine  that  is  based  on  the  desired  range  of  the  bins 

. */ 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  <math.h> 

#include  <fcntl.h> 

#include  <string.h> 

#define  DEBUG  1 
float  *VECT,  ‘HVECT; 

FILE  *infile,  *outfile,  *binfile; 


void  *  getvect(int  n,int  size) 
{ 


void  ‘vectptr; 


vectptr  =  (void  *)  calloc(n,size); 
if(vectptr  ==  NULL)  { 
printfC'Error:  dynamic  memory  fail\n"); 
exit(O): 

} 

return  (vectptr); 

} 


void  intquant(float  *vector,float  *headvector,int  size, 
int  stddev,int  mean,int  debug,int  numbins) 


float  fstddev,  fmean,  sd3  ,binsize: 
int  i=0; 


headvectorjO]  =  (float)mean; 
headvector[li  =  (float)stddev; 

fstddev  =  ((float)stddev)/100.0; 
fmean  =  ((float)  mean)/100.0; 
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/*  put  most  in  char  range  */ 
sd3  =  fstddev*3; 
binsize=  sd3/numbins; 

while  (i<size) 

{ 

if  (abs(vector[i]-fmean)<binsize) 

{ 

vector[i]  =  0; 

} 

else  vector[i]  =  (vector[i]-fmean)/binsize; 
i++; 

} 

} 

void  intunquant(fioat  ‘vector, float  *headvector,int  size.int  debug, int  numblns) 

Int  i=0,  stddev,  mean=0: 

float  fstddev,  fmean,  temp,  sd3,  binsize,  halfbin; 

mean  =  (int)headvector[0]: 
stddev  =  (int)headvector[1]: 

fstddev  =  ((float)stddev)/100.0; 
fmean  =  ((float)mean)/100.0; 
sd3  =  fstddev*3; 
binsize=  sd3/numbins; 
halfbin  =  binsize/2; 

printfC  mean  =  %6.2f,  stddev  =  %6.2f,  3sd  =  %6.2f,  binsize  =  %6.2f\n", 
fmean,fstddev,sd3,binsize); 

while  (i<size) 

{ 

if(vector[l]!=0)  vector[l]=(vector[i]*binslze)+fmean; 

else  vector[i]=fmean; 

i++; 

} 


void  quant{float  *vector,float  *headvector,int  size, 
int  stddev,int  mean, int  debug, int  numbins) 

float  fstddev,  fmean,  frtsdv,sd3; 

Int  i=0; 

fstddev  =  ((float)stddev)/100.0: 
fmean  =  ((float)  mean)/100.0; 
sd3  =  fstddev*3: 
frtsdv  =  sd3/numbins; 


headvector[0]  =  (float)mean; 
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headvector[1]  =  (float)stddev; 

while  (i<size) 

{ 

if  (abs(vector[i]-fmean)<frtsdv) 

{ 

vector[i]  =  0; 

} 

else  vector[i]  =  (vector[i]-fmean)/frtsdv; 
i++; 

} 

} 


void  unquant(float  ‘vector, float  *headvector,int  size,int  debug, int  numbins) 

{ 

int  i=0,  stddev,  mean=0,  rebuilt; 
float  fstddev,  fmean,  temp; 
float  sd3,  binsize,  halfbin; 


mean  =  (int)headvector[0]; 
stddev  =  (int)headvector[1]; 
fstddev  =  ((float)stddev)/100.0; 
fmean  =  ((float)mean)/100.0; 
sd3  =  fstddev*3; 
binsize  =  sd3/numbins; 
halfbin  =  binsize/2; 


printfC'mean  =  %6.2f,  stddev  =  %6.2f,  3sd  =  %6.2f,  binsize  =  %6.2f\n", 
fmean,fstddev,sd3,binsize); 

while  (i<size) 

{. 

if(vector[i]!=0)  vector[i]=(vector[i]*binsize)+fmean; 

else  vector[i]=fmean; 

i++; 

} 

} 


void  sectquant(float  *vector,float  ‘headvector,  int  size,int  zeroout,int  quantize,int  debug,int 
intqnt.int  numbins) 

{ 

int  i; 

double  high=-1 0000000.0,  low=1 0000000.0; 
int  zeros=0,  zrun=0,  zcount=0,  outsiders=0; 

double  sum=0.0,  squares=0.0,  mean=0.0,  var=0.0,  stddev=0.0,  sd3=0.0; 
double  rtsdv=0.0; 
float  binsize; 


if  (quantize==-1) 

{ 

if(intqnt)  intunquant(vector,headvector, size, debug, numbins); 
else  unquant(vector,headvector,size,debug, numbins); 
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} 

for  (i=0:i<size:i++) 

{ 

high  =  ((double)vector[i]>high)  ?  (double)vector[i] :  high; 
low  =  ((double)vector[i]<low )  ?  (double)vector[i] :  low  ; 
sum  +=  (double)vector[i]; 

squares  +=  (double)vector[i]*(double)vector[i]/(double)size; 

} 

mean  =  sum/(double)size; 
var  =  squares  -  (mean*mean); 
stddev  =  sqrt(var); 
sd3  =  3*stddev; 

binsize  =  sd3/numbins;  /*  zero  threshold,  3*std  dev  /bins  */ 

for  (i=0:i<size:i++) 

{ 

if  ((double)abs(vector[i]-mean)>sd3)  outsiders++; 
if  ((double)abs(vector[i]-mean)<binsize) 

{ 

zeros++; 

} 

else  if(zeroout) 

{ 

zeros++; 

vector[i]=0.0; 

} 

if  ((double)abs(vector[i]-mean)<binsize) 

{ 

zcount++; 

} 

else  zcount  =  0; 

zrun  =  (zcount>zrun)  ?  zcount :  zrun; 

} 

if  (quantize==1 )  printf ("  mean  =  %4. 1  f,  var  =  %8.2f,  stddev  =  %6.2f\n'', 
mean,  var,  stddev): 

printfC  high  =  %6.2f,  low  =  %6.2f,  size  =  %8i,  bins  =  %i  \n", 
high,  low,size,numbins); 

printf("  binsize  =  %6.2f,  zeros  =  %i,  longest  zero  run  =  %i\n", 
binsize,zeros,zrun); 

printf{"  outsiderthresh  =  %6.2f,  outsiders  =  %i\n",sd3, outsiders): 

if  (quantize==T) 

{. 

if  (intqnt)  intquant(vector,headvector,size,(int)(stddev*100), 
(int)(mean*100),debug,numbins): 
else  quant(vector,headvector,size,(int)(stddev‘1 00), 
(int)(mean*100),debug,numbins): 

} 

} 


void  main(int  argc,  char  *argv[]) 

{ 


132 


int  i,  quantize,  quant,  zero=0; 

int  xdim=0,  ydim=0,  size=0,  level=0,  sect=0,  cursize=0; 

int  numbins; 

char  clevel; 

float  f  level; 

float  *vector; 

float  ‘headvector; 

char  string[25]={0,0,0,0,0,0, 0,0,0, 0,0,0, 0,0,0, 0,0, 0,0,0,0,0,0,01,  *suf=''.bq'': 

if(argc  <  4) 

{ 

printfC'Usage  pgm  <filename>  Q  {1/-1}  \n"): 
printf(”First  name  is  source  file\n“); 
printf("Q=Levels  to  quantize  \n"); 

printf(''1  quantize,  -1  unquantize,  0  none\n''): 
exit(O); 

} 

infile  =  fopen(argv[1],"rb''); 
if(infile  ==  NULL) 

{ 

printfC'File  open  failure\n''); 
exit(O); 

} 

strcpy(string,argv[1]): 
strcat(string,".q2''); 
outfile  =  fopen(string,"wb"); 
if{outfile  ==  NULL) 

{ 

printfC'File  open  failure\n"); 
exit(O); 

} 

level  =  atoi(argv[2]); 
quantize  =  atoi(argv[3]); 
if  (abs(quantize)>1)  quantize  =  1; 
size  =  512*512; 

vector  =  (float  *)getvect(size,sizeof(float)); 

if  (quantize==1) 

{ 

fread(vector,sizeof(float),size,infile); 
headvector=(float  *)getvect(2*(level+1  ),sizeof(float)); 
flevel  =  (float)level; 

} 

else 

{ 

fread(&flevel,sizeof(float),1, infile); 
printf("%i  levels  to  process\n",(int)flevel); 
level  =  (int)flevel; 

headvector=(float  *)getvect(2*(level+1  ),sizeof(float)); 
fread(headvector,sizeof(float),2*(level+1), infile); 
fread(vector,sizeof(float),size,infiie); 

} 


for  (i=level;i>0;i--) 
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{ 

cursize  =  size/(int)pow(2.0,(double)i); 
if  (level==i) 

{ 

numbins=100; 

printfC'For  level  %i,  sect  %i  int  quantizing,  %i  bins  \n",  level, 1,numbins); 
sectquant(vector,&headvector[(level-i)*2],cursize,0,quantize,1,1,numbins); 

} 

for  (sect=2;sect<=2;sect++) 

{ 

numbins=(int)(100.0/flevel*(float)i); 

printfC'For  level  %i,  sect  %i  quantizing,  %i  bins\n",i,sect,numbins); 
sectquant(&vector[cursize],&headvector[(level-i+1)*2], 
cursize,0,quantize,1,0,numbins): 

} 

} 

if  (quantize==1) 

{ 

fwrite(&flevel,sizeof(float),1,outfile); 

fwrite(headvector,sizeof(float),(level+1)*2,outfile): 

fwrite(vector,sizeof(float),size,outfile); 

free(vector): 

free(headvector); 

} 

else 

{ 

free(headvector): 

fwrite(vector,sizeof(float),size,outfile); 

free(vector); 

} 

fclose(infile); 

fclose{outfile): 

} 
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. V 

FILENAME:  FSTATS.C 

PURPOSE:  Calculates  the  distribution  of  an  input  char  file 

Input  file  is  a  char  file 
output  file  is  a  text  file 


#include  <stdio.h> 
#include  <math.h> 
#include  <string.h> 
#include  <stdlib.h> 

#define  MAXLINE  512 
#define  NR_END  1 
#define  FREE_ARG  char* 
#define  MC  512 
#define  MQ  (2*MC-1) 


long  filesize(FILE  *stream) 

{ 

long  curpos,  length; 
curpos  =  ftell(stream): 
fseek(stream,  OL,  SEEK_END); 
length  =  ftell(stream): 
fseek(stream,  curpos,  SEEK_SET); 
return  length; 

} 


int  getfilename(char  name[25],int  argc,char  *argv[]) 

{ 

FILE  *in; 

char  destination[25],filename[25]; 

char  string[25]; 

int  i,namelength,error=1 ; 

while(error) 

{ 

for  (i=0;i<25;i++) 

{ 

string[i]=NULL; 

filename[i]=NULL; 

name[i]=NULL; 

} 

if(argc==1) 

{ 


printf("usage:  prgm  filename "); 
exit(O): 

} 

else 

{ 

strcpy(string,argv[1]); 
argc=1 ; 

} 

printfC'opening  %s\n", string); 
if  ((in  =  fopen(string,  "rb"))  ==  NULL) 

{ 

fprintf(stderr, "Cannot  open  %s  file,  press  return  \r", filename); 
error=getchar(); 

fprintf(stderr,"  \r"); 

error=1 ; 

} 

else 

{ 

fclose(in); 

error=0; 

} 

} 

return  1 ; 

} 


*ie*-k******ititifk*it***itifk-kit 


^  -  STATISTICS  ACCUMULATION  DRIVER 

int  main(int  argc.char  *argv[]) 

{ 

long  nfreq[256]; 

FILE  ’out,  ’in,  ’temp,  ’swap; 

FILE  ’tbi; 
char  filename[25], 

name[25HO,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}, 
ext[4]={0,0,0,0}. 
headerbuf[200]; 
unsigned  int  i,j,  getsize; 
long  fsize=0,runlen=0,count,charcount; 
unsigned  long  ’table=NULL; 
double  numbits; 
int  index; 

unsigned  char  inchar,  last; 
int  compress=0; 

float  entropy=0.0,prob=0.0,inf_num=0.0,inf_sum=0.0,avg_len=0.0,  sigma2=0.0; 

float  log2=log(2); 

float  mean=0.0,  meansq=0.0; 

double  sqsum=0.0,maxenergy=0.0; 


for(i=0;i<25;i++)  filename[i]=NULL; 
strcpy(filename,argv[1]): 

printfC'opening  in  %s\n'', filename); 

if  ((in  =  fopen(filename,"rb"))  ==  NULL) 

{ 

fprintf(stderr, "Input  %s  not  foundAn", filename); 
exit(1); 

} 

strncat(filename,".sts",4); 

printfC'opening  out  %s\n'', filename); 

if  ((tbi  =  fopen(filename,"w"))  ==  NULL) 

{ 

fprintf(stderr, "Output  file  %s  not  createdAn", filename); 
exit(1); 

} 

printf(“Analyzing  file  and  building  statistics\n"); 

fsize=filesize(in); 

for  (j=0;  j<256;  j++) 

{ 

nfreq[j]=0; 

} 

count=0; 

charcount=0; 

last=255;  /*  put  last  out  of  range  */ 

for(;;)  /*  character  counts*/ 

{ 

inchar=(unsigned  char)fgetc(in); 
index=((int)inchar);  /*  +128;*/ 
if(feof(in))  break; 
nfreq[index]++; 
count++; 

} 


fprintf(tbl,"ind  nfreq  prob  inf  num  \n\n"); 
printfC'ind  nfreq  prob  inf  num  \n\n"); 
for  (j— 0;  j<25G;j++) 

{ 

if  (nfreqlj]) 

{ 

prob=  (float)nfreq[J]/count; 
inf_num=  -(Iog(prob)/log2); 
inf_sum+=  inf_num; 
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sqsum  +=  j*j*nfreqO]; 
entropy+=  prob*inf_num: 
mean+=(float)j*prob; 

printf(''%3i,  %6i,  %6.4f,  %6.4f  \n",j,nfreq[j],prob,inf_num): 
fprintf(tbl,"%3i,  %6i,  %6.4f,  %6.4f  \n",j,nfreqO],prob,inf_num); 

} 

} 


meansq  =  mean’mean; 

for  (j=0:  j<256;j++) 

{ 

sigma2+=  (float)nfreq[j]*(((float)j-mean)*((float)j-mean)): 

} 

sigma2  =  sigma2/(float)(count); 

fprintf(tbl,"\n\nFILE_STATISTICS\n"); 
fprintf(tbl,"total_char=  %lu\n'', count); 
fprintf(tbl,"info_number=  %f\n'',inf_sum); 
fprintf(tbl,"entropy=%f\n'', entropy): 
fprintf(tbl,"mean=  %An",  mean); 
fprintf(tbl,"variance=  %An'',sigma2): 
fprintf(tbl,"std_dev=  %An'',sqrt((double)sigma2)); 
fprintf{tbl,"totaLenergy=  %An“, sqsum); 

printfC'ENCODING  STATISTICS\n"): 
printf("total  char  read  =  %lu\n", count); 

printf("information  number  =  %An'',inf_sum); 
printf("file  entropy 
printfC'mean 
printf(”variance 
printf(“std  deviation 
printf("total_energy 

fclose(tbl); 
fclose(in); 
return  0; 

} 


=  %An", entropy): 

=  %An",  mean); 

=  %An",sigma2); 

=  %An",sqrt((double)sigma2)): 
=  %An", sqsum); 


/* . . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  FCMP.C 

PURPOSE:  Builds  a  comparison -error  distribution 

Input  files  are  char  files 
output  file  is  a  text  file 


#include  <stdio.h> 
#include  <math.h> 
#include  <string.h> 
#include  <stdlib.h> 

#define  MAXLINE  512 
#define  NR_END  1 
#define  FREE_ARG  char* 
#define  MC  512 
#define  MO  (2*MC-1) 


long  filesize(FILE  *stream) 

{ 

long  curpos,  length; 
curpos  =  ftell(stream): 
fseek(stream,  OL,  SEEK_END); 
length  =  ftell(stream); 
fseek(stream,  curpos,  SEEK_SET); 
return  length; 

} 


MAIN  -  STATISTICS  ACCUMULATION  DRIVER 


★  *****************^**'*****************i^**************imr******************** 


int  maln(int  argc,char  *argv[]) 

{ 

long  nfreq[256]; 

FILE  *out,  *in,*in2,  *temp,  ‘swap; 

FILE  *tbl; 
char  filename[25], 

name[25]={0,0,0,0,0,0,0,0,0,0, 0,0,0, 0, 0,0,0,0,0,0,0,0,0,0,01, 
ext[4]={0,0,0,0}, 
headerbuf[200]; 
unsigned  int  i,j,  getsize; 
long  fsize=0,runlen=0,count,charcount; 
unsigned  long  *table=NULL; 
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double  numbits; 
int  index; 

unsigned  char  inchar,  inchar2,  last; 
int  compress=0; 

float  entropy=0.0,prob=0.0,inf_num=0.0,inf_sum=0.0,avg_len=0.0,  sigma2=0.0; 

float  log2=log(2); 

float  mean=0.0,  meansq=0.0; 

double  sqsum=0.0,maxenergy=0.0; 

for(i=0;i<25:i++)  filename[i]=NULL; 
strcpy(f  ilename,argv[1  ]); 

printfC'opening  in  %s\n", filename); 

if  ((in  =  fopen(filename,''rb"))  ==  NULL) 

{ 

fprintf(stderr, "Input  %s  not  found.\n",&argv[1][0]); 
exit(1); 

} 

if  ((in2  =  fopen(argv[2],"rb"))  ==  NULL) 

{ 

fprintf(stderr, "Input  %s  notfound.\n”,&argv[2][0]); 
exit(1); 

} 

strncat(filename,".dif'',4); 

printfC'opening  out  %s\n", filename); 

if  ((tbi  =  fopen(filename,"w"))  ==  NULL) 

{ 

fprintf(stderr, "Output  file  %s  not  createdAn", filename); 
exit(1); 

} 

printfC’Analyzing  file  and  building  statisticsVn"); 

fsize=filesize(in); 

for  (j=0;  j<25S;  ]++) 

{ 

nfreq[j]=0; 

} 

count=0; 

charcount=0; 

last=255;  /*  put  last  out  of  range  7 

for(;;)  /*  character  counts7 

{ 

inchar=(unsigned  char)fgetc(in); 
inchar2=(unsigned  char)fgetc(in2); 
index=((int)inchar-inchar2+128);  /*  +128;7 
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if(feof(in))  break; 

nfreq[index]++: 

count++: 

} 


fprintf(tbl,"ind  nfreq  prob  inf  num  \n\n"): 
printf("ind  nfreq  prob  inf  num  \n\n''); 
for  (j=0:  j<256:j++) 

{ 

if  (nfreqUl) 

{ 

prob=  (float)nfreq[j]/count; 

inf_num=  -(log(prob)/log2): 

inf_sum+=  inf_num; 

sqsum  +=  0-128)*(j-128)*nfreq0]: 

entropy+=  prob*inf_num; 

mean+=(float)(j)*prob; 

printf(''%3i,  %6i,  %6.4f,  %6.4f  \n",j-128,nfreq|]],prob,inf_num); 
fprintf(tbl,''%3i,  %6i,  %6.4f,  %6.4f  \n",j-128,nfreq|j],prob,inf_num); 

} 

meansq  =  mean*mean; 
for  (j=0;  j<256;j++) 

{ 

^  sigma2+=  {float)nfreq[j]*(((float)j-mean)*((float)j-mean)); 

sigma2  =  sigma2/(float)(count): 

fprintf(tbl,"\n\nFILE_COMPARISON_STATISTICS\n"); 
fprintf(tbl,''%s  -  %s\n\n",&argv[1][0],&argv[2][0]): 
fprintf(tbl,''totaLchar=  %lu\n", count); 
fprintf(tbl,''info_number=  %f\n",inf_sum); 
fprintf(tbl,"entropy=  %f\n", entropy); 
fprintf(tbl,''mean=  %f\n'',  mean-128.0); 
fprintf(tbl,"variance=  %f\n",sigma2); 
fprintf(tbl,”std_dev=7of\n",sqrt((double)sigma2)); 
fprintf(tbl,''totaLenergy=  7of\n'', sqsum); 


printfC'DIFFERENCE  STATISTICSNn”); 
printf("%s  -  %s\n\n'',argv[1],argv[2]); 
printf("total  char  read  =  7olu\n", count); 

printf("information  number  =  7of\n",inf_sum); 
printf("file  entropy  =  7of\n", entropy); 


printf("mean 
printf("variance 
printf("std  deviation 
printf("total_energy 


=  7of\n",  mean-128.0); 

=  7of\n",sigma2); 

=  %f\n",sqrt((double)sigma2)); 
=  %An", sqsum); 


fclose(tbl); 
fclose(in); 
fclose(in2); 
return  0; 

} 


/* . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  RGBLOG.C 

PURPOSE:  Splits  a  24-bit  RGB  file  into  three  planes  using  a 
transform  based  on  the  human  visual  system 

Input  file  is  a  24-bit  color  file 
3  char  output  files  created  containing  a  plane  each 

NOTE:  similar  to  RGBKKK  and  KKKRGB  see  notes  for  inverse  transform 


. . */ 

#include  <stdio.h> 

#include  <math.h> 

#include  <string.h> 

#include  <stdlib.h> 

r  file  size  determinator  */ 
long  filesize(FILE  *stream) 

{ 

long  curpos,  length; 
curpos  =  ftell(stream): 
fseek(stream,  OL,  SEEK_END): 
length  =  ftell(stream): 
fseek(stream,  curpos,  SEEK_SET): 
return  length; 

} 


^**********^**-*******|^^|^****************************^*****^**y 

int  main(int  argc,char  *argv[]) 

{ 

FILE  *in,  ‘R,  *G,  *B,  *Y,  *1,  *Q,  *YT,  *IT,  *QT; 
int  i,j,bmp=0; 

char  filename[1 2],othername[1 2]={0,0,0,0,0,0,0,0,0,0,0,0}, 
name[8]={0, 0,0,0, 0,0,0,0}, 
ext[4]={0,0,0,0}, 
headerbuf[100]; 

unsigned  char  RGB[3],KKK2[3]; 
float  KKK[3]; 
int  width,height; 
unsigned  long  count=0; 
long  size; 
int  textout=0; 

if  (argc<2)  printf("usage:  pgm  filename  [t]  {for  text  output}\n"); 

strcpy(name,argv[1]); 

if  ((in  =  fopen(argv[1],"rb"))  ==  NULL) 


{ 

fprintf(stderr, "Input  %s  not  foundAn", filename); 
exit(1); 

} 

size=filesize(in): 

height=width=(int)sqrt((double)size/3): 

if((long)height*width*3-size!=0) 

{ 

printfC'filesize  error\n“); 
exit(1): 

} 

printf("height=%i,  width=%i,  size=%lu\n", height, width, (unsigned  long)height*width); 

for  (i=0;i<=14;i++)  othername[i]=0; 
strcpy(othername,name); 
strncat(othername,".l1  ”,4); 

if  ((Y  =  fopen(othername,"wb”))  ==  NULL) 

{ 

fprintf(stderr,"Output  file  %s  not  createdAn",othername); 
exit(1); 

} 

for  (i=0;i<=14:i++)  othername[i]=0; 

strcpy(othername,name); 

strncat(othername,".l2",4): 

if  ((I  =  fopen(othername,"wb"))  ==  NULL) 

{ 

fprintf(stderr, "Output  file  %s  not  createdAn",othername); 
exit(1): 

} 

for  (i=0;i<=14:i++)  othername[i]=0; 

strcpy(othername,name); 

strncat(othername,".l3",4); 

if  ((Q  =  fopen(othername,"wb"))  ==  NULL) 

{ 

fprintf(stderr,''Output  file  %s  not  createdAn",othername); 
exit(1); 

} 

for(i=0;i<height*width;i++) 

{ 

fread(RGB,1,3,in): 

KKK[0]=(0.299*(float)RGB[0]+0.587*(float)RGB[1]+0.114*(float)RGB[2]); 

KKK[1]=(0.607*(float)RGB[0]+0.174*(float)RGB[1]+0.20r(float)RGB[2]); 

KKK[2]=(  0.066*(float)RGB[1]+1.117*(float)RGB[2]); 


KKK[0]=log(KKK[0]): 
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KKK[1]=log(KKK[1]): 

KKK[2]=log(KKK[2]); 

KKK2[0]=(unsigned  char)((21.50*KKK[0]+106.277)/.6197); 
KKK2[1]=(unsigned  char)((-41 .0*KKK[0]+41 .00*KKK[1]+21 .652)/.1344): 
KKK2[2]=(unsigned  char)((-6.27‘KKK[0]+6.27*KKK[2]+1 7.771  )/.0941); 


fwrite(&KKK2[0],1,1,Y); 
fwrite(&KKK2[1], 1,1,1); 
fwrite(&KKK2[2],1,1,Q): 

} 

printf("\n"); 

fclose(in): 
fclose(Y); 
fclose(l): 
fclose(Q); 
return  0; 

} 
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/* . . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  VF.C 


PURPOSE;  Float  file  viewer  for  the  SGI  machines,  displays 
values  from  0  -  255  truncates  values  beyond  that  to  0  or  255 

Input  file  is  a  binary  float  file 

NOTE:  8-bit  greyscale  filesare  similarly  handled  just  make  the 
appropriate  changes  to  handle  and  display  char  based  files 


#include  <gl/gl.h> 
#include  <gl/device.h> 
#include  <stdio.h> 
#include  <fcntl.h> 
#include  <math.h> 


7 


void  *  getvect(int  n,int  siz) 

{ 

void  *ptr; 


ptr  =  (void  *)  calloc(n,siz): 
if(ptr  ==  NULL)  { 

printf("Error:  Dynamic  memory  allocation  fail\n"); 
exit(O): 

} 

return(ptr); 

} 


/*  file  size  determinator  7 
int  fiiesize(FILE  ‘stream) 
{ 


int  curpos,  length; 
curpos  =  ftell(stream): 
fseek(stream,  OL,  SEEK_END); 
length  =  ftell(stream): 
fseek(stream,  curpos,  SEEK_SET); 
return  length; 

} 


main(int  argc,char  ‘argvQ) 

{ 


intx; 
int  y; 
FILE  *fp; 
int  i; 
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int  n; 
int  X=0: 
int  Y=0; 
int  C; 

int  waiting  =1 ; 
int  notdone  =1 ; 

Device  dev; 
int  size,hsize=0; 
float  scaler=1 .0; 
int  offset=0; 
short  colo[256][3]; 
long  vert[2]; 
float  *c; 

float  ‘image,  ‘header; 
short  val; 

if(argc  <  2)  { 

printfC'This  program  is  for  float  format  files  \n"): 
printfC'Usage:  pgm  filename  [options]\n"); 
printf("-s#  scaler  =  multipler  for  values  float#\n“); 
printf(''-o#  offset  =  offset  values  by  int#\n“); 
printf("-h#  size  of  header  \n"); 
printf("-x#  force  width  to  int#\n"); 
printf("-y#  force  width  to  int#\n"): 
printf("-b  windows  bmp  format  {not  implemented}\n"); 

printf("if  the  file  is  one  of  the  standard  sizes  it  will  automatically  determine  the  size\n"); 
printfC'in  program:  up  &  down  arrow  keys  scale  ‘2r.5\n"); 
printff  left  &  right  arrow  keys  offset  -20l+20\n"): 
exit(O); 

} 

fp=fopen(argv[1],''rb''): 
if  (fp==NULL) 

{ 

printfC'file  open  error\n"): 
exit(O); 

} 

for  {i=2;i<argc;i++) 

( 

switch  (argv[i][1]) 

{ 

case  ■x':X  =  atoi(&argv[i][2]); 
break; 

case  'y':Y  =  atoi(&argv[i][2]); 
break; 

case  'o':offset  =  atoi(&argv[i][2]); 
break; 

case  's'lscaler  =  atof(&argv[i][2]); 
break; 

case  'h':hsize  =  atoi(&argv[i][2]); 

header  =  (float  *)  getvect(hsize  ,sizeof(float)); 

fread(header,sizeof(float),hsize,fp); 

break; 
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} 

} 

size  =  filesize(fp)/si2eof(float): 

switch  (size-hsize) 

{ 

case  262144:X=512; 

Y=512: 

break; 

case  65536  ;X=256; 

Y=256; 

break; 

case  64000  :X=320; 

Y=200; 

break; 

case  16386  :X=128; 

Y=128; 

break; 

case  4096  :X=64 ; 

Y=64; 

break; 

default  :if(X>0&&Y==0) 

Y=(size-hsize)/X; 
else  if(Y>0&&X==0) 
X=(size-hsize)A'; 
else  if(Y+X==0) 

X=Y=(int)sqrt((double)size-hsize); 

} 

image  =  (float  *)  getvect(X*Y,sizeof (float)); 

prefsize(X,Y); 

winopen(argv[1]); 

ortho2(0.0,X-1,-Y+1,0.0); 

cmodeO; 

gconfigO; 

color(WHITE); 

clearO; 

qdevice(ESCKEY); 

qdevice(REDRAW); 

qdevice(RIGHTARROWKEY); 

qdevice(LEFTARROWKEY); 

qdevice(UPARROWKEY); 

qdevice(DOWNARROWKEY); 

n=fread(image,sizeof(float),X*Y,fp); 

for(x=1  ;x<256;x++)  { 
mapcolor(x,x,x,x); 

} 

gconfigO; 

while(notdone) 

{ 
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c=image: 
for(y=0;y<Y;y++)  { 
for{x=0:x<X;x++)  { 

if  ((int)(*c+offset)*scaler>255)  C  =  255; 
else  if  ((int)(*c+offset)*scaler<  0)  C  =  0; 
else  C  =  (int)(*c+offset)*scaler: 
color(C): 

C++; 

bgnpointO; 

vert[0]=x; 

vert[1]=-y; 

v2i(vert); 

endpointO; 

gflushO; 

} 

} 

waiting=1 ; 
while(waiting) 

{ 

while  (!  qtest()){} 
switch(dev=qread(&val)) 

{ 

case  REDRAW:  waiting=0; 
break; 

case  UPARROWKEY: 
case  DOWNARROWKEY: 
case  RIGHTARROWKEY: 
case  LEFTARROWKEY: 
case  ESCKEY:  if(val==0)  waiting=0; 
break; 

} 

} 

switch{dev) 

{ 

case  ESCKEY  :  notdone  =  0;  break; 
case  UPARROWKEY: 
scaler*=2: 

printf  ("scaler  set  to  %5.2f\n", scaler); 
break; 

case  DOWNARROWKEY: 
scaler/=2; 

printf  ("scaler  set  to  %5.2f\n", scaler); 
break; 

case  RIGHTARROWKEY:  offset  +=  20; 

printf  ("offset  set  to  %3i\n'', offset); 
break; 

case  LEFTARROWKEY:  offset  -=  20; 

printf  ("offset  set  to  %3i\n", offset); 
break; 

case  REDRAW : 


} 

/*  while(! (qtestO  &&  qread(&val)  ==  ESCKEY  &&  val  ==  0))  { }7 


} 

free(image); 
free(header); 
gconfigO: 
gexitO: 
return  0; 

} 


/* . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  VCOLOR.C 

PURPOSE:  Displays  24-bit  rgb  images  on  the  SGI  machines 
Input  file  is  a  24-bit  rgb  file  also  will  handle  windows  BMP  files 


#include  <gl/gl.h> 
#include  <gl/device.h> 
#include  <stdio.h> 
#include  <fcntl.h> 
#include  <math.h> 


*/ 


void  *  getvect(int  n,int  siz) 
{ 


void  *ptr; 

ptr  =  (void  *)  calloc(n,siz); 
if(ptr  ==  NULL)  { 

printf("Error:  Dynamic  memory  allocation  fail\n“); 
exit(O); 

} 

return  (ptr): 

} 


★********/ 


/*  file  size  determinator  V 
int  filesize(FILE  ‘stream) 

{ 

int  curpos,  length; 
curpos  =  ftell(stream); 
fseek(stream,  OL,  SEEK_END): 
length  =  ftell(stream); 
fseek(stream,  curpos,  SEEK_SET): 
return  length; 

} 


main(int  argc.char  *argv[]) 

{ 
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int  x; 
int  y: 

FILE  ‘fp: 
int  i,j; 
int  n; 
int  X=0: 
int  Y=0: 
float  C; 
int  bmp=0; 
int  waiting  =1 ; 
int  notdone  =1 ; 

Device  dev; 
int  size,hsize=0; 
float  scaler=1.0: 
float  offset=0; 
short  colo[256][3]; 
long  vert[2]; 

float  blackvect[3]  =  {0.0, 0.0, 0.0}; 
float  colovect[3]; 
unsigned  char  *c; 
unsigned  char  *image,  ‘header; 
short  val; 

if(argc  <  2)  { 

printf("This  program  is  for  unsigned  char  format  files  \n"); 
printfC'Usage:  pgm  filename  [options]\n“); 
printf("-s#  scaler  =  multipler  for  values  float#\n“); 
printf("-o#  offset  =  offset  values  by  float#\n"); 
printf("-h#  size  of  header  \n“); 

printf("-x#  force  width  to  int#\n''); 
printf("-y#  force  width  to  int#\n"); 
printf(”-b  windows  bmp  format  \n''); 

printf("if  file  is  one  of  the  standard  sizes  it  will  automatically  determine  the  sizeXn"); 
printfC'in  program:  up  &  down  arrow  keys  scale  *2r.5\n''); 
printf("  left  &  right  arrow  keys  offset  -.05l+.05\n"); 
exit(O); 

} 

fp=fopen(argv[1],"rb”); 
if  (fp==NULL) 

{ 

printfC'file  open  errorXn"); 
exit(O); 

} 

size  =  filesize(fp)/(sizeof(unsigned  char)); 

for  (i=2;i<argc;i++) 

{ 

switch  (argv(i][1]) 

{ 

case  'x':X  =  atoi(&argv[i][2]); 
break; 

case  y:Y  =  atoi(&argv[i][2]); 


break; 

case  'o':offset  =  atof(&argv[i][2]); 
break; 

case  's':scaler  =  atof(&argv[i][2]); 
break; 

case  'h'lhsize  =  atoi(&argv[i][2]); 

header  =  (unsigned  char  *)  getvect(hsize  ,sizeof(unsigned  char)); 

fread(header,sizeof(unsigned  char),hsize,fp); 

break; 

case  'b':hsize  =  54; 

header  =  (unsigned  char  *)  getvect(hsize  ,sizeof(unsigned  char)); 
bmp  =  1; 

fread(header,sizeof(unsigned  char),hsize,fp); 
break; 

} 

} 

size  -=  hsize; 
size  /=  3; 


switch  (size) 

{ 

case  262144:X=512; 

Y=512; 

break; 

case  65536  :X=256; 

Y=256; 

break; 

case  64000  :X=320; 

Y=200; 

break; 

case  16386  :X=128; 

Y=128; 

break; 

case  4096  :X=64 ; 

Y=64: 

break; 

default  :if(X>0&&Y==0) 

Y=(size)/X; 
else  if(Y>0&&X==0) 

X=(size)A'; 
else  if(Y+X==0) 

X=Y=(int)sqrt((double)size); 

} 

printfC'image  is  %i  bytes  with  a  %i  header,  X=%i  by  Y=%i\n", size, hsize, X,Y); 

image  =  (unsigned  char  *)  getvect(size*3,sizeof(unsigned  char)); 

prefsize(X,Y); 

winopen(argv[1]); 

ortho2(0.0,X-1,-Y+1,0.0); 

RGBmodeQ; 

gconfigO; 

color(WHITE): 

c3f(blackvect); 


clear(): 

qdevice(ESCKEY); 

qdevice(REDRAW); 

qdevice(RIGHTARROWKEY): 

qdevice(LEFTARROWKEY): 

qdevice(UPARROWKEY): 

qdevice(DOWNARROWKEY); 

n=fread(image,sizeof(unsigned  char),size*3,fp); 

gconfigO: 

while(notdone) 

{ 

c=image; 

for(y=0;y<Y;y++) 

{ 

for(x=0;x<X;x++) 

{ 

if  (bmp) 

{ 

for  (j=2:j>=0:j-) 

{ 

if  ((((float)*c/255.0)*scaler)+offset>=1.0)  C  =  0.99; 
else  if  ((((float)*c/255.0)*scaler)+offset<  0.0)  C  =  0.0; 
else  C  =  (((float)*c/255.0)*scaler)+offset; 
colovect|j]=C; 

C++; 

} 

c3f(colovect); 

bgnpointO; 

vert[0]=X-x; 

vert[lj=-Y+y; 

} 

else 

{ 

for  O=0;j<3;j++) 

{ 

if  ((((float)*c/255.0)*scaler)+offset>=1 .0)  C  =  0.99; 
else  if  ((((float)*c/255.0)*scaler)+offset<  0.0)  C  =  0.0; 
else  C  =  (((float)*c/255.0)*scaler)+offset; 
colovect[j]=C; 

C++J 
}  • 

c3f(colovect); 

bgnpointO: 

vert[0]=x; 

vert[1]=-y: 

} 

v2i(vert): 

endpointO; 

gflushO; 

} 


waiting=1 ; 
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while(waiting) 

{ 

while  (!  qtest()){} 
switch(dev=qread(&val)) 

{ 

case  REDRAW:  waiting=0: 
break; 

case  UPARROWKEY: 
case  DOWNARROWKEY: 
case  RIGHTARROWKEY: 
case  LEFTARROWKEY: 
case  ESCKEY:  if(val==0)  waiting=0: 
break; 

} 

} 

switch(dev) 

{ 

case  ESCKEY :  notdone  =  0;  break; 
case  UPARROWKEY: 
scaler*=2.0; 

print!  ("scaler  set  to  %5.2f\n“, scaler); 
break; 

case  DOWNARROWKEY: 
scaler/=2.0; 

printf  ("scaler  set  to  %5.2f\n'', scaler); 
break; 

case  RIGHTARROWKEY:  offset  +=  .05; 

printf  ("offset  set  to  %3.0f\n'', offset); 
break; 

case  LEFTARROWKEY:  offset  -=  .05; 

printf  ("offset  set  to  %3.0f\n", offset); 
break; 

case  REDRAW  : 


} 

/*  while(l (qtestO  &&  qread(&val)  ==  ESCKEY  &&  val  ==  0))  { }7 

} 

free(image); 
free(header); 
gconfigO; 
gexitO; 
return  0; 

} 
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/* . 

THOMAS  W.  PIKE,  Thesis  Project,  Spring  '95 


FILENAME:  VIEWRGB.C 
PURPOSE:  PC  based  color  image  viewer 


■V 


*  THOMAS  W.  PIKE,  SPRING  95,  VERS  2.0 

* 


*  * 

PURPOSE:  RGBcolor  graphics  viewer  for  SVGA  systems  * 

*  * 

FEATURES:  automatic  file  size  detection  routines,  builds  a  256  psuedo 

*  color  map  for  24  bit  images 

*  * 

*.  NOTES:  compiled  with  turbo  C++  source  option  on  * 

*  * 

*********************************^***************^^*******,j.************* ! 

#include  <stdlib.h> 

#include  <graphics.h> 

#include  <stdio.h> 

#include  <conio.h> 

#include  <math.h> 

#include  <svga256.h> 

#include  <string.h> 

#define  MAX_WIDTH  2048 


/*  */ 

/*  FUNCTION  FILESIZE  returns  the  filesize 

/*  7 

/*  */ 


V 


long  filesize(FILE  'stream) 
{ 


long  curpos,  length; 
curpos  =  ftell(stream); 
fseek(stream,  OL,  SEEK_END): 
length  =  ftell(stream); 
fseek{stream,  curpos,  SEEK_SET): 
return  length; 

} 


/*  7 

/*  PROCEDURE  GETFILENAME 
/* 

/* 


7 

7 


7 


int  getfilename(char  ‘name,  char  *ext,  int  argc.char  *argv[]) 

{ 

FILE  ‘in; 

char  destination[25],f  ilename[1 2]; 
char  string[20],  ‘dot; 
int  i,namelength,error=1; 

while(error) 

{ 

for  (i=0;i<20;i++) 

{ 

string[i]=NULL: 
if(i<12)  filename[i]=NULL; 
if(i<8)name[i]=NULL: 
if(i<4)  ext[i]=NULL; 

} 

if(argc==1)  { 

printfC'viewrgb  filename  [-f  -r]"); 
printf(''option  f  flips  the  order  of  display"); 
printf("option  r  flips  the  rgb  order  to  bgr"); 


} 

else 

{ 

strcpy(string,argv[1]): 
argc=1 ; 

} 

if  ((in  =  fopen(string,  "rb"))  ==  NULL) 

{ 

fprintf(stderr, "Cannot  open  %s  file,  press  return  \r", filename); 
error=getchar(); 

fprintf(stderr,"  \r"); 

error=1 ; 

} 

else 

{ 

fclose(in); 

error=0; 

} 

} 

dot=strchr(string,'.'); 
namelength  =  dot-string; 
strncpy(name,  string,  namelength); 
for(i=namelength;i<8;i++)  name[i]=NULL; 
strncpy(ext,  dot,  4); 


return  1; 
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r 


*/ 


/*  DetectVGA256  is  a  dummy  function  for  the  TURBO  C  initialization 


int  huge  DetectVGA256(void) 

{ 

return  0; 

} 


/*  */ 

/*  FUNCTION  setsvga256  graphics  initialization 
/*  7 


/* 


7 


void  setsvga256{void) 
/* - 


GRAPHICS  INITIALIZATION  BLOCK 


- */ 

{ 

/*  select  the  svga256  driver  and  mode  that  supports  the  use  7 
/*  of  the  setrgbpalette  function.  7 

int  gdriver,  gmode,  errorcode; 
int  ht,  xmax; 

/*  Install  the  new  driver  7 

installuserdriver("SVGA256",DetectVGA256): 

gdriver=DETECT; 

/*  initialize  graphics  and  local  variables  7 
initgraph(&gdriver,  &gmode, 

/*  read  result  of  initialization  7 
errorcode  =  graphresult(): 

if  (errorcode  !=  grOk) 

/*  an  error  occurred  7 

{ 

printfC'Graphics  error:  %s\n",  grapherrormsg(errorcode)); 

printf("Press  any  key  to  halt:”); 

getchO: 

exit(1);  /*  terminate  with  an  error  code  7 

} 

/* - 


END  OF  GRAPHICS  INITIALIZATION  BLOCK 
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^**********************************Hf<f****************************<f  ******  j 

/*  */ 

/*  FUNCTION  SETColorSCALE  creates  256  psuedo  colorscale 
/*  7 

/*  7 

^***********************************************imr**<r*****************^i(lf  j 

void  setcolorscale(void) 

{ 

int  rcolor,gcolor,bcolor,ncolor: 
int  i=0,j=0,k; 

/*  create  color  scale  7 

for  (rcolor=0;  rcolor<8;  rcolor++) 

{ 

for  (gcolor=0;  gcolor<8;  gcolor++) 

{ 

for  (bcolor=0:  bcolor<8;  bcolor++) 

{ 

ncolor=((rcolor«6)+(gcolor«3)+bcolor): 

if(ncolor%2==0) 

{ 

ncolor/=2; 

printf("%7i\r“,ncolor): 

setrgbpalette(ncolor,rcolor*8,gcolor*8,bcolor*8); 

} 

} 

} 

} 

r  display  palette  7 
for(i=0;i<256:i++) 

{ 

for(j=0:j<8;j++) 

{ 

putpixel(i*2,520+j,i): 
putpixel(r2+1 ,520+j,i): 

} 

} 

} 

^*****************************************it,t*************************y 

void  main(int  argc,  char  *argv[]) 

{ 

/*  initialization  7 
long  fsize; 
int  i,j,k,l,  ask=1 ; 

FILE  *  infile; 
int  A.B.C,  color; 


7 


char  buffer[MAX_WIDTH],  inchar; 
unsigned  char  COLOR,  OFFSET; 
char  filename[1 2],name[8],ext[4]: 
int  header,  width,  height,  offset; 
float  scaler; 
int  flip=0,  reverse=0; 

/*  open  file  V 

if(!getfilename(name,ext,argc,argv)) 

{ 

fprintf(stderr,"lnput  file  not  found,  terminating  programV); 
exit(O); 

} 

if  (argc>2) 

{ 

for(i=2;i<argc;i++) 

{ 

if(argv[i][0]='-') 

{ 

switch(argv[i][1]) 

{ 

case  'f:  flip=1 ; 
break; 

case  'r':  reverse=1 ; 
break; 

} 

} 

} 

} 

for  (i=0;i<12;i++)  filename[i]=NULL; 
for  (i=0;i<8;i++)  filename[i]=name[i]; 
strncat(filename,ext,4); 
if  ((infile  =  fopen(filename,"rb''))  ==  NULL) 

{ 

fprintf(stderr,"lnput  %s  not  found.\n'',filename); 
exit(1); 

} 

fsize=filesize(infile)/3; 

/*  calculate  file  parameters  */ 
if(fsize>=64000&&fsize<66048) 

{ 

width=320; 
height=200; 
header=fsize-64000; 
if  (header==0)  ask=0; 

} 

else  if (fsize>=1 28000&&fsize<1 28200) 

{ 

width=640; 
height=200; 
header=fsize-1 28000; 
if  (header==0)  ask=0; 


} 

else  if (f s  ize>=224000&&f s  ize<224200) 

{ 

width=640; 
height=350; 
header=fsize-224000; 
if  (header==0)  ask=0; 

} 

else  if(fsize>=:245760&&fsize<245960) 

{ 

width=512: 
height=480; 
header=fsize-245760: 
if  (header==0)  ask=0: 

} 

else  if{fsize>=307200&&fsize<307400) 

{ 

width=800: 
height=600; 
header=fsize-307200; 
if  (header==0)  ask=0; 

} 

else  if(fsize>=745560&&fsize<745760) 

{ 

width=(int)(sqrt((double)fsize)); 

height=width; 

header=fsize-width*height; 

} 

else 

{ 

width=(int)(sqrt((double)fsize)); 

height=width; 

header=fsize-(width*height); 

header*=3: 

} 

if  (ask) 

{ 

printfC'Header  size  =  %i\n'', header): 
printfC'lmage  width  =  %i\n'', width); 
printfC'lmage  height  =  %i\n”, height); 
printf("Are  these  correct?  [y.n]"); 
scant  ("%c",&lnchar) : 
printf("\n"): 

} 

if  (inchar=='n') 

{ 

printfC'lmage  width  = "); 

scanf(“%i’',&width); 

printfC'lmage  height  = "); 

scanf("%i",&height): 

header=fsize-height*width; 

header*=3: 
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printf("Header  estimate  =  %i\n'', header); 
printf("Header  size  = "); 
scanf(“%i“,&header); 

} 

if  (argc  >=  3)  offset=  (int)  atol(argv[2]); 
else  offset=0; 

if  (argc  >=  4)  scaler=  atof(argv[3]); 
else  scaler=1 .0; 

/* . 


GRAPHICS  INITIALIZATION  BLOCK 


- */ 

setsvga256(); 

setgraphmode(3);  /*  mode  3  is  800X600  */ 
setcolorscaleO: 

/* . 


END  OF  GRAPHICS  INITIALIZATION  BLOCK 


. 7 

/*  file  display  routine  7 

fread(&buffer,  sizeof(char),  header,  infile); 

for  (i=  0;  i<height;  ++i) 

{. 

if  (fread(&buffer,  sizeof(char),  width*3,  infile)==width*3) 

{ 

for  0=0;  j<width*3;  j+=3) 

{ 

r  scaler  multiplies  range  by  value  to  spread  out  or  contract  range 
offset  changes  the  starting  point  of  the  color  range,  anything 
pushed  above  255  with  the  shift  will  be  left  at  255 
7 

if(reverse)  COLOR  =(((buffer[j+2]/32)«6)+((buffer[j+1  ]/32)«3)+(buffer[j]/32))/2; 
else  COLOR  =(((buffer[j]/32)«6)+((buffer[j+1  ]/32)«3)+(buffer0+2]/32))/2; 
if(flip)  putpixel(width-j/3,height-i,(COLOR)); 
else  putpixelO/3,i,(COLOR)); 

} 

} 

else  printfC'full  buffer  not  read"); 

} 

r  wait  until  done  viewing  7 
getchO; 

/*  cleanup  and  exit  7 
fclose(infile); 
closegraphQ; 
restorecrtmodeO; 

} 
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